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1 . INTRODUCTION 


A theory  of  designing  the  "best"  processor  for  estimating  the  state  of 
a Markov  process  or  signal  observed  indirectly  as  a zero  memory  nonlinear 
function  of  the  state,  corrupted  by  additive  white  noise,  was  developed  in  the 
early  and  middle  1960s.  "Best"  here  is  to  mean  the  estimate  that  minimizes  the 
expected  losfl:  with  the  loss  a function  of  the  estimation  error.  This  theory 
is  a generalization  of  the  Kalman-Bucy  linear  filtering  theory.  However,  unlike 
the  Kalman-Bucy  theory,  no  blueprint  for  the  black  box  which  accepts  as  inputs 
the  observations  and  produces  as  outputs  the  optimal  estimate,  the  minimum  loss 
estimate,  can  be  deduced  from  the  theory.  In  fact,  in  general,  the  black  box, 
nonlinear  filter  corresponding  to  the  optimal  estimate  is  infinite  dimensional 
in  that  the  state  of  the  nonlinear  filter  is  not  finite  dimensional.  For 
important  technological  problems,  the  signal  process  estimators  generated  by 
linearization  and  application  of  the  linear  Kalman-Bucy  design  often  exhibit 
poor  performance,  sometimes  even  divergence.  Further,  even  when  the  linear 
design  seems  effective,  one  is  interested  in  what  estimator  has  the  best 
possible  performance  and  what  this  performance  is,  in  order  to  know  whether  further 
effort  on  this  optimal  design  is  justified.  Finally,  study  of  the  optimal 
estimator  allows  one  to  generate  effective  suboptimal  designs. 

For  the  above  reasons,  it  became  clear  that  building  nonlinear  filters 
was  an  important  task.  Now  it  is  quite  easy  to  see,  for  example  see  [10], 
that  the  black  box,  determining  the  optimal  filter,  is  specified  by  the  maps 
with  domain  at  time  n,  the  observation  process  sample  path  up  to  and 
including  time  n and  range,  the  conditional  probability  of  the  signal  at 
time  n given  the  observation  process  sample  path  up  to  and  including  n, 

Jn(.,  z^, zj,  with  Zj  the  observations. 


Our  problem  then  has  two  phases,  one  the  replacement  of  Jn(*.  ^ z^) 

by  a finite  vector  J^,  the  representation  problem,  and  second,  the  realiza- 
tion problem  which  is  the  means  used  to  generate  J^i  from  and  z +, . 

The  realization  tool  which  has  been  used  has  been  digital  and  hybrid 
systems.  The  problems  considered  have  been  the  cubic  sensor,  passive  receiver 
and  the  one,  two  and  three  state  dimensional  phase  demodulation  problem. 

These  problems  have  been  chosen  for  their  importance,  but  also  because  the  state 
dimension  was  low,  a necessity  if  sufficiently  precise  estimates  of  the  error 
performance  of  the  optimal  filter  are  to  be  generated.  This  is  because 
error  performance  can  so  far  only  be  evaluated  by  Monte  Carlo  simulation. 

This  is  also  true  for  suboptimal  filters.  The  problems  which  we  have  studied 
the  most  are  the  cubic  sensor  problem  and  the  two-state  dimensional  phase 
demodulation  problem;  see  [4],  [11],  [12],  & [13]. 

The  representations  considered  for  the  two-dimensional  phase  demodulation 
problem  have  been:  point  mass,  Gauss-Hermite  polynomials,  Fourier  Series,  and 
Cubic  Splines  under  tension. 

On  the  basis  of  the  structure  of  the  equations  to  be  solved  in  order 
to  build  an  optimal  nonlinear  filter,  it  became  clear  early  in  our  synthesis 
research  effort  that  parallel  or  associative  digital  computers  offered 
considerable  speed-up  in  estimate  production  over  standard  serial  machines 
of  third  generation;  see  [1]  and  [4]  where  these  effects  are  discussed. 

Because  of  our  lack  of  access  to  machines  of  the  parallel  or  associative 
type  in  the  early  1970's,  a parallel  processor  was  constructed  using  a 
contemporary  hybrid  system.  In  [2]  the  feasibility  was  considered,  while 
in  [3]  , [5]  the  construction  of  the  actual  system  was  reported  along  with 


subsequent  Monte  Carlo  runs,  all  for  a one-dimensional  signal  process  cubic 
sensor  problem.  The  results  were  encouraging;  in  particular  the  hybrid 
realization  proved  to  be  16  times  faster  than  an  all-digital  realization, 
using  the  digital  portion  of  the  hybrid  system.  Our  ultimate  aim,  however,  was  to 
study  the  effect  of  using  CDC  Star  100  and  Illiac  IV  as  synthesis  tools  for 
the  nonlinear  filter.  This  paper  is  devoted  to  detailing  the  impact  of  these 
machines  on  the  nonlinear  filtering  problem  of  phase  demodulation. 

In  the  summer  of  1975,  we  obtained  access  to  the  Illiac  IV  and  made 
plans  to  access  the  CDC  Star  100,  when  it  was  delivered  to  NASA  Langley 
Research  Center.  It  was  decided  that  in  order  to  effectively  be  able  to 
judge  the  improvement,  we  should  initially  realize  a problem,  the  two-dimen- 
sional phase  demodulation  problem  where  we  had  extensive  numerical  results 
and  experience  on  the  serial  CDC  6600.  Further  for  this  problem  we  could 
generate  estimate  and  signal  sequences  with  an  existing  serial  6600  program. 

We  were  convinced  that  this  serial  program  was  close  to  the  fastest 
possible,  and  that  that  was  not  the  case,  was  an  interesting  byproduct  of 
producing  an  effective  Star  program.  While  studying  this  two-dimensional 
phase  demodulation  problem  on  the  Illiac  and  the  Star  gave  us  good  data  on 
speed-up  factors  possible  with  pipeline  and  parallel  machines,  our  real 
objective  was  to  do  a problem  which  was  beyond  the  capabilities  of  third 
generation  serial  machines.  Consequently,  we  developed  a generalized  three- 
state  dimension  demodulation  problem  where  phase,  phase  rate  and  amplitude 
all  must  be  estimates.  In  the  Fall  and  the  Spring  of  1975-1976,  Ken  Senne 
developed  a two-dimension  optimal  phase  demodulator  program  for  the  Illiac  IV 
using  the  Glypnir  language  and  the  Arpa  net.  In  the  month  of  June,  the  authors 
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were  visiting  scientists  at  ICASE,  NASA  Langley.  During  this  time,  the  corre- 
sponding program  for  the  Star  100  was  developed.  This  latter  program  not  only 
was  very  effective  for  Star  (a  large  percentage  of  operations  were  of  long 
vector  type,  over  4,000  component  vectors,  consequently  assuring  that  Star 
was  doing  long  streaming),  but  also  this  Star  program  was  used  as  a prototype 
of  an  extremely  fast  serial  program  which  was  over  2-1/3  times  faster  on  the 
6600  than  our  original  program  reported  in  [4];  also  see  the  listing  of  the 
cyclic  point  mass  program.  Finally,  a Star  program  for  a three-dimensional 
signal  process,  combining  amplitude  and  phase  estimation,  was  developed  and  time 
comparison  between  it  and  the  corresponding  serial  version  obtained. 

Methods  of  obtaining  extremely  fast  serial  realization  of  the  optimal 
filter  are  being  investigated  using  optical  and  surface  wave  devices.  It  is 
curious  that  the  mathematical  mapping  which  permits  a time  correlator  surface 
wave  device  to  calculate  multi-dimensional  convolutions  is  intrinsic  in  the 
Star  program— compare  [9].  We  will  discuss  later  in  this  paper  the  possible 
application  of  array  processors  coupled  with  minicomputers  to  this  problem. 

In  [10],  a discussion  is  given  comparing  the  speed-up  possibilities 
inherent  in  realization  method  versus  density  representation.  This  paper  is 
concerned  only  with  realization  methods  and  uses  essentially  only  the  point 
mass  representation.  Other  representations  can  be  used,  but  in  general  they 
must  be  accuracy  calibrated  against  the  point  mass  method, and  the  effect  of 
pipeline  and  parallel  realization  on  estimate  speed-up  is  unclear. 


2.  COMPUTER  REALIZATION  SURVEY 


2.1  The  Physical  Problem 

It  has  been  observed  in  previous  studies  ([4])  that  nonlinear  filters 
can  be  efficiently  implemented  on  parallel  computers.  Since  every  so-called 

i 

vector  machine  has  unique  limitations,  however,  it  Is  necessary  to  adapt 
the  algorithm  to  each  particular  architecture.  Candidate  architectures  for 
the  present  study  include  the  array  processor  (Illiac  IV),  the  pipeline 
processor  (COC-Star  100),  and  the  look-ahead  processors  (CDC  6600  and  7600). 
Benchmarks  on  fast  serial  machines  (IBM  370-168  and  POP  11-70)  have  also  been 
included  for  reference. 

An  interesting  application  for  optimal  nonlinear  estimation  was 
introduced  by  Malllnckrodt,  Bucy,  and  Cheng  [7],  who  considered  the  problem 
of  tracking  a first-order  phase  process  based  on  measurements  of  a modulated 
signal  in  noise  of  the  form 

ds(t)  * A cos  [u)0t  + x^tjjdt  + dv(t)  (2.1) 

where  A is  a known  amplitude,  wQ  is  a known  carrier  frequency,  and 
xj(t)  is  the  message  process  being  tracked.  The  measurement  noise  is 
assumed  white.  Using  a voltage-controlled  oscillator,  the  known  carrier 
may  be  removed  by  heterodyning  down  to  base  band,  producing  both  in-line 
and  quadrature  components  and  resulting  in  an  equivalent  two-dimensional 
measurement  process  of  the  form 


'57 


where  A has  been  taken  as  unity  without  loss  of  generality,  and  the 
noise  has  been  replaced  by  a vector  of  mutually  independent  quantities. 

The  first-order  phase  process  studied  in  [7]  consisted  of  Brownian 
motion  with  increment  of  length  h having  variance  qh.  In  this  paper, 
we  describe  a study  of  a second  order  phase  process  involving  the  integral 
of  Brownian  motion,  expressed  as 


0 1 X] 

0 0 Xo 


(2.3) 


We  will  retain  the  same  measurement  model  (1)  and  let  the  noises  v^  and 
v^  be  independent  Brownian  motions  with  paths  of  length  h having  variance 
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The  optimal  filter  is  specified  by  Jn+i i n+i  , the  conditional 

i y 2 

density  of  the  phase  and  phase  rate  at  time  (n  + 1)a  , given  discrete 
observations  up  to  this  time  at  sampling  interval  a.  The  following 


equations  determine  this  density: 


j fc  W,4 


(yp-w)  ) /y-j -yA 

-Mm  , 


° co exp  { 


z-j(n+l)  cos  y-j  + z2(n+l)  sin  y1 
r/A 


(2.5) 


It  can  be  shown- -see  [4]—  that  a modulated  density  Jn|n  defined  as 
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with  -n<_a  < it,  and  -ir/A  <_  t < n/A,  carries  all  the  information  necessary 
for  nonlinear  filtering  when  the  cyclic  loss  function  ^(1  - cos  e-j ) with 
the  phase  error,  is  used.  The  modulated  density  satisfies 


where 


We  recall  that  the  cyclic  estimate,  the  one  that  minimizes  the  cyclic  loss 


2.2  Point  Mass  Placement 


The  most  attractive  formulation  for  numerical  solution  of  the  filtering 
problem  on  parallel  machines  involves  the  use  of  point-mass  representation  of 
the  densities.  A fixed  rectangular  grid  will  consist  of  m points  in  the 
phase  variable  and  n points  in  the  phase  rate,  with  coordinates  (i,  j) 
corresponding  to  (a,s)  by  the  formulae 


Note  that  the  phase  rate  variables  may  be  used  directly  in  the  convolution 


[2  1)  but  the  phase  variables  must  be  interpolated  in  general,  so  that 


n(k)  = o(i)  - A?  (j)  =>  k = i + | - J j + j) 


(2.12) 


If  m is  taken  as  an  even  integer,  such  that  n/m  is  an  integer,  then 

I 

(2.12)  may  be  decomposed  into  a subdominant  integer  [k]  and  a fractional 
part  Ak,  as  follows:* 


[k(i )]  = 1 +|  - ( j + j ) DIV  (£)  - 1 . 


(2.13) 


where  DIV  denotes  integer  division.  Note  that  the  addition  of  — in  (2.13) 
will  not  change  the  result  since  n/m  is  an  integer.  Next,  we  observe 
that  the  remainder  after  the  division  in  (2.13)  is  the  fractional  part  Ak; 


( 1 + j)  MOD  ( - ),  j - 0,  ....  n-1 
2 m 


(2.14) 


Ak  - 1 - ( j + j ) MOD  ( “■  ) . 


The  interpolated  result  desired  will  be  between  [ k ( i ) ] and  [ k ( i ) ] +1 , 
(evaluated  modulo  m)  with  weighting  Ak  on  the  former  and  1 - Ak  on  the 
latter.  If  the  convolution  is  formed  for  fixed  phase  (i),  then 


[k(i  + 1)]  = [k(i ) ] + 1 (module  m)  . 


(2.15) 


The  relation  (2.15)  suggests  the  recursion  which  will  be  used  on  the  Illiac. 
(2.15)  will  be  used  to  iterate,  based  on  an  initial  condition 


[k(  )]  = ( - - 1 - j DIV  “ ) MOD  m 

i v /j  '2  J m 


Note  that  [k]  should  be  interpreted  as  modulo  m. 


(2.16) 
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During  such  iteration,  i,  Ak  will  be  used  to  weight  the  term  evaluated 
at  [ k ( i ) ] , while  1 - Ak  will  be  used  to  weight  the  term  evaluated  at 
[k(l  + 1)]. 
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5(j) 
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(2.17) 


2.3  Evaluating  the  Filtered  Data  Term 

The  term  (2.8)  consists  of  an  infinite  sum.  The  values  which  r 
may  take  on  in  the  grid  coordinates  consist  of  integer  multiples  of  2ir/nA, 
resulting  in 

a(r(p) ) - l exp  f-  ~ ( j)  ( n * i )J  p - -2(n-l).  — ,2(n-I) 
t='"  (2.18) 

It  may  be  seen  from  (2.18)  that  a(*)  is  an  even  function  which  is  cyclic, 
modulo  n.  Further,  for  reasonable  values  of  q the  contribution  of  all  terms 
except  £ = 0,  is  negligible.  Thus,  we  compute 


a(r(p))  = exp 
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qA 


m (M) 
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(2.19) 
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In  the  examples  described  below,  5 terms  of  a(*)  were  taken  to  be  non-zero 
for  the  computing  convolutions. 


2.4  Assessment  of  Required  Computations 


The  computations  required  to  implement  the  point  mass  filter  are 
summarized  in  Table  2.1,  as  a function  of  m and  n.  The  sensor  terms 
are  the  only  ones  which  require  math  functions  (exponentials).  Since  only 
m exponentials  are  required,  this  computation  is  generally  insignificant 
compared  with  the  overall  filter  update,  so  no  special  effort  has  been 
expended  to  optimize  the  required  computations. 


TABLE  2.1  ARITHMETIC  OPERATIONS  FOR  FILTER 


Function 

Multiplies 

Number  of  Operations 
Adds  Divisions 

Exponentials 

Sensor  Terms 

2m 

m 

0 

m 

Interpolation 

mn 

2mn 

0 

0 

Convol ution 

5mn 

lOmn 

0 

0 

Row  Sums 

0 

mn  m 

0 

0 

Estimates 

3m 

3m- 3 

1 

0 

Normalization 

mn+m 

0 

0 

0 

Total 

7mn+6m 

13mn+3m-3 

1 

m 

Example 

m = 32 
n a128 

28864 

53341 

1 

32 

2.5  The  Illiac  IV  Algorithm 


The  primary  considerations  for  programming  the  Illiac  IV  array  are 
the  proper  utilization  of  all  of  the  64  Processor  Elements  (PE's)  and 
minimization  of  data  routing  between  PE's.  In  order  to  accomplish  the  efficient 
use  of  the  PE's,  the  values  of  m = 32  and  n = 128  were  selected  and  utilized 
for  all  machine  examples.  On  the  Illiac,  two  rows  of  PE  storage  were  used 
for  each  value  of  the  phase  samples  in  JN.  The  interpolation  and  convolution 
were  accomplished  in  a totally  parallel  fashion* using  all  64  PE's,  except 
for  the  operation  depicted  in  Figure  2.1.  It  is  necessary  to  perform  a 
cyclic  rotation  of  each  phase  row  to  accumulate  the  terms  for  the  convolution 
(a  cyclic  rotation  to  the  right  is  combined  with  a cyclic  rotation  to  the 
left  at  each  step).  Since  a cyclic  routing  on  Illiac  IV  involves  only  one 
row  at  a time,  it  was  necessary  to  form  the  two-row  rotation  by  two  single-row 
rotations,  followed  by  an  end  element  switch  (involving  three  transfers 
with  only  one  PE  enabled). 

The  overall  effectiveness  of  the  Illiac  IV  algorithm  is  evident  from 
the  fact  that  so  few  operations  are  required  which  involve  fewer  than  64  PE's 
enabled.  The  sensor  terms  are  computed  with  only  32  PE's  enabled,  but  this 
operation  only  involves  about  5%  of  the  estimate  update.  The  single  PE  trans- 
fers in  the  convolution  are  also  only  responsible  for  about  5%  of  the 
computation  time.  Finally,  the  row  sums  are  done  logarithmically  with  a 
PE  utilization  efficiency  of  16.7%,  and  they  account  for  another  5%  of  the 
estimate  computation  time.  Thus,  the  net  efficiency  of  this  algorithm  is 
87.5%,  or  the  equivalent  of  56  times  faster  than  a single  PE  program. 


The  Illiac  program  was  coded  in  the  Glypnir  language,  utilizing 
assembly  language  listings  to  reduce  the  overhead  computations.  The  resulting 


program  was  timed  in  the  non-overlap  mode,  wherein  the  Control  Unit  (CU) 
and  the  PE's  operate  serially  and  considerable  overhead  cycles  are  required 
to  allow  error  conditions  to  ring  out  between  instructions.  Thus,  although 
we  would  expect  the  Illiac  program  to  run  at  least  twice  as  fast  as  the 
equivalent  program  for  the  CDC-STAR,  in  fact  it  ran  five  times  slower.  It 
will  be  reported  in  [17]  how  the  program  operates  in  overlap  mode. 

It  is  also  important  to  indicate  that  the  Illiac  IV  Clock  is  running  at 
80  nsec,  as  compared  with  the  design  goal  of  50  nsec. 

2.6  The  CDC-STAR  Algorithm 

The  pipeline  architecture  is  unconstrainted  by  small  fixed  resources 
(i.e.,  64  processors).  On  the  other  hand,  efficient  utilization  of  the 
pipeline  requires  detailed  attention  to  pre-arrangement  of  vectors  to 
allow  for  streaming  from  consecutive  memory  locations.  This  consideration  is 
particularly  important  for  STAR,  which  has  a relatively  slow  memory  cycle 
time.  Since  the  nonlinear  filter  is  recursive,  it  is  necessary  to  include 
the  vector  re-arrangement  as  part  of  the  filter  update  and  therefore  the 
re-arrangement  constitutes  the  major  overhead  of  the  STAR  program.  The 
operations  on  the  matrix  JN  to  precondition  it  for  the  convolution  are  shown  in 
Figures  2.2  and  2.3.  First,  the  column-ordered  JN  matrix  is  column-shuffled 
with  itself  to  produce  a matrix  which  has  two  copies  of  every  phase  variable  in 
each  column.  Then,  a scrambled  JN  can  be  formed  which  has  the  property 
that  each  row  in  the  final  convolved  matrix  can  be  generated  by  operating 
on  a suitable  interpolation  between  the  two  adjacent  rows  of  the  scrambled  JN. 
The  interpolation  which  does  this  is  depicted  in  Figure  2.3. 


*The  Figures  2.3  - 2.4  are  shown  with  m = 4 and  n = 16  for  illustrative 
purposes  only. 


Figure  2.3.-  Interpolation  of  scrambled  matrix. 


This  interpolation  may  be  done  by  vector  operations  of  length  (m+l)n,  or  4224. 
The  row  corresponding  to  m+1  in  the  result  may  be  compressed  out  of  the  final 
result  to  reduce  the  subsequent  calculations. 

The  cyclic  convolution  is  shown  in  Figure  2.4.  First,  the  end  columns 
of  the  interpolated  JN  are  cyclically  copied  to  produce  a matrix  from  which 
each  of  the  terms  of  the  convolution  may  be  obtained  as  mxn  matrices. 

The  production  of  a 5- term  symmetric  convolution  is  done  in  parallel 
for  all  4096  points  by  a sequence  of  10  vector  adds  and  5 vector  multiplies, 
all  of  length  4096. 

The  only  computations  which  are  done  on  the  STAR  that  are  less  than 
100%  efficient  are  the  vector  sums  and  the  determination  of  the  estimates. 

This  is  reflected  in  Table  2.2,  which  gives  the  breakdown  of  the  various 
functions  in  the  STAR  program. 

2.7  CPC  6600  Program 

The  CDC  6600  (and  7600)  has  instrurtion  look-ahead  and  multiple  arithmeti 
processor  units  which  provide  a partial  overlap  parallelism.  The  most 
efficient  6600  programs  contain  many  tight  loops  instead  of  complicated 
computations  within  large  loops.  By  recoding  the  vectorized  STAR  program 
in  analogous  FORTRAN  for  the  CDC  6600,  we  were  able  to  achieve  near  optimal 
utilization  of  the  available  resources  with  functional  loops  which  are  for 
the  most  part  contained  within  the  8-word  instruction  stack. 

The  basic  data  flow  of  the  CDC  6600  is  illustrated  in  Figure  2.5. 

Reads  from  memory  are  accomplished  by  setting  the  A Registers  A1  - A5  with 


Convolved  density. 
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TABLE  2.2  --  STAR  PROGRAM  BREAKDOWN 


Operations  % of  Total  start-ups  Compute  Required 


Vector  Arithmetic 

63.4 

5174 

76793 

78.5% 

12  Multiplies 

33.3 

1908 

41152 

70.1% 

46  Adds 

26.4 

3266 

30861 

86.4% 

1 Exponential 

3.7 

- 

4780 

100.0% 

Vector  Rearrangement 

31.3 

1233 

39292 

83.7% 

1 Vector  Transfer 

5.8 

1001 

6464 

68.0% 

2 Indexed  Transfers 

22.1 

144 

28480 

100.0% 

(Block  Lengths  32  & 33) 

1 Vector  Compress 

3.4 

88 

4348 

0% 

Scalar  Arithmetic 

79 

100.0% 

1 Divide 

46 

3 Adds 

33 

Subroutine  Overhead 

4.1 

5142 

0% 

3 Calls 

. _ . _ _ _ 

Miscellaneous 

1.2 

1537 

100.0% 

• 

Memory  Conflicts,  etc. 

- 

— 

TOTAL 

6407 

122843 

=>  5.17msec 

5% 

95% 

"■i  Minimum  Achievable 

6224 

111296 

=>  4.70msec 

the  appropriate  address;  writes  are  obtained  by  loading  A6-A7.  The  B 
Registers  are  used  for  incrementing  and  address  computation. 

The  breakdown  of  the  computations  for  the  CDC  6600  is  shown  in 
Table  2.3.  Note  that  the  major  overhead  is  for  reading  and  writing  and 
miscellaneous  waiting.  It  is  interesting  to  note  also  that  the  ratio  of 
multiply  rates  between  STAR  and  6600  is  25  to  1 , while  the  add  rate  ratio 
is  17.5  (including  normalization  on  6600).  Thus,  for  arithmetic  alone, 

STAR  would  be  expected  to  be  21  times  faster  than  the  6600  on  this  problem. 
The  achieved  speed-up  of  29  is  explained  by  the  slightly  lower  overhead  on 
STAR  (36.6%  versus  52.9%  for  6600). 


2.8  Implications  for  Future  Synthesis  of  Higher  Dimensional  Filters 


The  main  purpose  of  this  study  was  to  measure  the  speed-up  achievable 
by  parallel  and  pipeline  machines,  in  order  to  determine  the  feasibility  of 
realizing  3 and  4 state  dimensional  nonlinear  filters.  In  a later  section 
of  this  paper  we  describe  the  fairly  natural  extension  of  our  2-state 
dimensional  filter  STAR  100  program  to  a 3-state  dimensional  problem  of 
estimating  phase,  phase  rate  and  amplitude.  One  outgrowth  of  the  software 
development  was  a determination  of  machine  structure  limitations  for  natural 
synthesis  of  higher  dimensional  filters.  Basically,  we  found  that  Illiac  IV 
was  limited  to  the  two-dimensional  problem  because  of  limited  P.E.  memory 
and  the  availability  of  only  64  PE's,  and  most  importantly,  because  of  the 
increase  of  program  complexity  and  overhead  with  dimension.  The  STAR  100 
program  also  becomes  complex  to  program  and  overhead  prone,  but  at  a somewhat 
higher  dimension  with  the  basic  limitation  that  the  total  number  of  grid  points  be 
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TABLE  2.3-  CDC  6600  PROGRAM  BREAKDOWN 


Function 
(%  Cycles) 

Operations 

Multiplies  Adds  Divisions 

Reads 

Wri tes , 

Exp. 

Notes 

Sensor  Update 
('0) 

2m 

m 

0 

3m 

m 

m 

Not  in  stack 
(extn.  refs.) 

Interpolation 

(-16%) 

mn 

2mn 

0 

4mn+3n 

mn 

0 

Outer  Loop 
not  in  stack 

Expansion 

('0) 

0 

0 

0 

10m 

10m 

0 

Instack 

Convolution 

(-74%) 

5mn 

lOmn 

0 

16mn+15 

6mn 

0 

Instack 

Row  Sums 

(-4%) 

0 

mn-m 

0 

mn+m 

m 

0 

Instack 

Estimates 

(~0) 

3m 

3m- 3 

1 

4m+l  1 

3 

0 

Instack 

Normalization 

(-6%) 

mn+m 

0 

0 

mn+2m 

mn 

0 

Instack 

Approx.  Totals 

7mn 

I3mn 

0 

22mn 

8mn 

0 

Approx.  Number 
of 

Minor  Cycles 

70mn 

91  mn 

0 

* 

66mn 

24mn* 

0 

Summary: 
Minor  Cycles 
(Approx. ) 

Arithmetic: 

Read/Write: 

Waiting: 

161  mu 
90mn 
91mn 

47.1% 

26.3% 

26.6% 

Measured: 

342mn 

100.0% 

(m  = 32,  n = 

1 28  for 

test  case) 

Measured 

Time  = 

140  msec/est. 

★ 

To  some  extent  these  cycles  are  concurrent,  but  the  total  of  read/write 


less  than  65K,  the  maximum  length  of  vector  for  Star  vector  operations, 
as  we  show  later,  a 3-state  problem  is  feasible  on  Star. 


The  implications  for  future  machine  design  suitable  for  our  problem 
are  a pipeline  machine  which  can  compute  with  indexed  vectors  and  with 
reduced  operation  cycle  time. 

3.  TIME  FACTORS  FOR  2-DIMENSIONAL  PHASE  DEMODULATION 

In  [4]  a complete  description  of  the  two-state  dimensional  phase 
demodulation  problem  is  given  as  well  as  numerous  numerical  results  and 
a listing  of  the  cyclic  nonlinear  filtering  program.  Because  we  wished  to 
compare  the  Illiac,  Star  and  6600  directly,  the  grid  of  32*128  was  chosen, 
i.e.,  32  subdivisions  in  phase  and  128  subdivisions  in  phase  rate,  to 
represent  JN  the  conditional  density  of  phase  and  phase  rate  given  the 
observations.  This  grid  was  natural  for  the  Illiac  as  it  has  64  parallel 

channels  so  that  the  integral  equation 

-ir/A 

Dn(x)  / 

Jn(x,y)  = J S(y-s)  J ,(x-CA,c)  ds  (3.1) 

k -n/A  n_l 

y € (-it/ a,  Tt/ a)  X e ( -TT  ,TT  ) 

with  the  values  of  Jn(x<-  ) could  be  found  in  two  passes;  the  first 

simultaneously  giving  { J ( x • , y.)}  and  the  second  giving  simultaneously 

0 j=l . . .64 

{Jn(x-i,  y-)}  and  then  repeated  for  each  i.  We  knew  before 

" 3 j=65. ...  128 

that  this  mesh  was  fine  enough  for  accuracy.  The  serial  cyclic  point  mass 
filter  with  interpolation  required  .700  seconds  per  estimate. 


The  Illiac  on  the  other  hand  produced  estimates  every  .0308  second.  The 
Illiac  was  not  run  at  full  capacity  for  our  runs  and  we  were  unable  to 
run  it  under  the  fastest  possible  conditions  because  other  jobs  had  tied  up  the 
machine  in  August,  so  that  our  timing  should  be  looked  at  as  an  upper  bound 
to  actual  performance  , clearly  far  away  from  the  predicted  theoretical 
performance.  It  is  unclear  whether  in  the  production  mode  Illiac  can 
achieve  anywhere  near  theoretical  performance,  i.e.,  64  times  faster  than 
CDC  6600  on  a full  parallel  job.  It  may  be  that  Illiac  hardware  cannot  be 
made  to  operate  at  its  theoretical  speeds. 

As  described  in  the  previous  section  for  the  Star  100  program,  the 
density  was  carried  as  a 32  xl28  = 4096  component  vector  and  the  convolution 
developed  as  a sum  of  scalar  times  translates  of  this  vector.  This  was 
done  in  order  to  take  advantage  of  the  pipeline  speed  of  Star.  At  first 
with  a direct  Star  Fortran  program,  a rate  of  .0118  second  per  estimate  was 
obtained.  By  examination  of  the  assembly  language  version  of  the  Star  program 
and  timing  each  element,  we  replaced  certain  Star  Fortran  calls  by  more 
efficient  assembly  language  routines,  e.g,  VGATHER  and  VSUM  were  replaced. 

After  this  was  done,  we  achieved  a rate  of  .0049  second  per  estimate. 

From  this  latter  number  it  became  clear  that  the  serial  6600  program  was  not 
efficient.  We  were  exceeding  the  factor  of  25  theoretical  speed  ratio  of 
Star  to  the  6600.  This  ratio  is  the  ratio  of  Star  multiplication  time  to 
6600  multiply  time.  A new  program,  Starrun,  was  written  which  was  the 
serial  version  of  the  Star  program.  Starrun  achieved  .178  second  per  estimate 
with  the  FTN,  OPT  = 2,  Level  410  Fortran  compiler.  The  FTN  4.5 
compiler  produces  significantly  different  and  faster  assembly  language 
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programs  than  the  earlier  optimizing  compilers.  By  elimination  of  some 
multiplications  and  arranging  things  so  that  everything  is  in  the  stack, 
the  Starrun  program  was  made  to  run  at  .14  second  per  estimate.  The 
following  table  gives  the  timing  of  our  programs  on  various  machines. 
With  considerable  effort  in  each  case  devoted  to  generating  fast  code: 


Machines 

Milliseconds  Per  Estimate 

Compiler 

CDC  6600 

137.0 

FTN  4.5  Level  410  0PT=2 

STAR  100 

4.7 

Assemby  Language  Optimization 

IBM  370-168 

130.0 

H Compiler 

ILL I AC  IV 

30.8 

Glypnir  Language,  Assembly 
Listing  Optimization 

PDP  11-70 

870.0 

FORTRAN  4 plus 

CDC  7600 

25 

FTN  4.5  Level  410  0PT=2 

In  order  to  get  a feel  for  the  relative  speeds  of  the  Star  and  the 
6600,  the  ratio  of  multiply  times,  when  the  pipe  on  the  Star  is  filled  is  25  to  1 , 
i.e.,  6600  takes  10  minor  cycles  at  100  nano  seconds  per  cycle,  while  a 
multiply  on  Star  takes  1 cycle  at  40  nanoseconds  per  cycle.  However,  the 
6600  has  several  look-ahead  features  as  well  as  a stack  which  can  speed  up 
the  processing.  The  FTN,  0PT=2  compiler  seeks  to  put  all  loops  in  the  stack 
and  to  fill  all  the  registers  so  that  the  ratio  we  have  found  28.57  to  1 is 
very  reasonable  for  a problem  such  as  ours  which  is  structured  so  that  vector 
operations  can  be  done  most  of  the  time  on  Star.  Of  course,  using  Star  for 
a problem  which  does  not  have  parallel  structure  would  be  wasteful  as  one 
would  do  as  well  or  better  on  the  7600. 
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4.  BENCHMARK  APPLICATION  FOR  DIFFERENT  APPROXIMATIONS 


^ i 


It  is  very  popular  in  engineering  studies  to  come  up  with  a compromise 
realization  which  approaches  the  accuracy  of  the  optimum  nonlinear  filter  with 
tremendous  amount  of  computation  savings.  Without  the  existence  of  an  accurate 
optimum  scheme,  which  serves  as  a benchmark  against  which  any  other  approxima- 
tion schemes  can  be  judged,  little  can  be  done.  The  development  of  such  approxi- 
mations is  possible  when  the  exact  performance  and  amount  of  computations  is 
known.  Then  and  only  then,  the  engineering  compromise  can  be  studied  intel- 
ligently. 

One  of  these  compromise  studies  has  recently  been  published--see  [14]-- 
as  an  approximation  for  the  optimum  two-dimensional  phase  demodulation.  It 
can  be  generalized  to  three-dimensional  phase-amplitude  demodulation  problems. 
This  approximation  is  based  on  fitting  the  probability  density  function  at 
each  time  step  to  a Gaussian  distribution  which  produces  the  optimum 
first  and  second  moments  for  the  same  error  criterion.  This  means  that  the 
update  is  accurate  and  required  only  for  the  first  and  second  moments. 

In  order  to  be  able  to  investigate  the  extension  of  this  Gaussian  approximation 
to  the  three-dimensional  problem,  it  is  essential  to  study  very  carefully  the 
behavior  of  the  optimum  filter  and  establish  a performance  benchmark. 

This  3-dimensional  problem  was  neither  analyzed  in  the  optimum  way  nor 
simulated  on  any  digital  computer  before  because  of  its  size  and  complexity. 

4 first  attempt  was  carried  out  to  realize  the  optimum  demodulator  on  the 
CYBER  175  and  CDC  STAR-100  in  order  to  establish  boundary  limits  on  the 
size  of  the  program  storage  and  the  computation  time. 

The  mathematical  model,  optimum  filter  and  preliminary  results  for  the 
3-dimensional  problem  will  be  discussed  in  the  following  secions. 


’’TT '"V""  
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5.  INTRODUCTION  TO  30-PROBLEM  AND  PRELIMINARY  RESULTS 


Achieving  a good  speed-up  factor  for  the  two-dimensional  phase 
demodulation  problem  was  the  signal  to  move  one  step  ahead  and  do  the  three- 
dimensional  problem  (phase,  phase  rate,  and  amplitude  demodulation). 


5.1  Mathematical  Model 

The  three-dimensional  state  vector  consists  of  phase  x,  phase  rate  y 
with  additive  white  Gaussian  noise  u and  amplitude  A corrupted  with  another 
independent  white  Gaussian  noise  v additive  to  the  amplitude  mean  value. 


Vi  + Vi  A 


Vi  + Vi 

a(a  - An-1 )a  + A„_1  + vn_1 


where  a is  the  time  step 

a is  the  damping  constant 
a is  the  amplitude  mean  value. 

E u =0 
n 

E = q,4 
£»„  = ° 

E V*  = q24 
E*0  ■ 0 

2 _ \ 3/** 
E x0  = / 2 q1  r 

e y0  = 0 
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The  observation  z is  modelled  in  the  following  way; 


l . l 

z = A cos  x + wn 
n n n n 


2 . 2 

z = A„  sin  xn  + w_ 
n n n n 


12 


1 2 

, 


wn  is  white  Gaussian  noise  independent  of  un,  vn  with  E wn  = 0 


1,2  2 

E(w_  ) 


r_ 

A 


A is  chosen  to  be  the  same  as  in  the  two-dimensional  problem,  i.e.. 


_ r % 

■'  [* ( ?,>  J 


c1  is  chosen  to  be  2.19  (same  as  2-dimensional  problem) 


c^  is  chosen  to  be  4 


3 1 


a = 1 


q2  - -01 


5.2  Optimum  Filter 


The  optimum  filter  realization  is  carried  out  by  a recursive  update 
of  the  density  Jn  using  the  discrete  form  of  the  presentation  theorem 


--see  [4]. 


f f'h 

Jn(x,y,A)  = Kn  Sn(x,A)  J J a^y-z)  a2(A-Br)-aaA)j  , (x-zA,z,n)dzdn 

— tt/ A 


2 ^ 

where  Sn(x,A)  = Exp  jj  [A(z„  cos  x + z*  sin  x)  - | ]J 

a-j (y  -z)  - l Expf--1-  (y  - z ♦ ^ )'] 

1 k=-«  L 2q-| A J 

a?(A  -6n-aaA)  = Exp f-  — — (A  -en-aaA)  1 

L 2q2A  J 


S = 1 - aaA 


K is  normalizing  constant  chosen  so  that 
n 


f f [ jn(x,  y.  A)  d x dy  dA  = 1 
-®  — tt/ A -ir 
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5.3  Optimum  Estimates 


The  phase  estimates  are  produced  in  the  same  manner  as  in  two- 
dimensional  problems.  They  are  obtained  by  minimizing  the  error  function 
L](e)  = 2(1  - cos  e) 


(x,  y.  A)  dx  dy  dA 


The  amplitude  estimates  are  produced  by  simply  minimizing  the  quadratic 

o 

loss  function  l^e)  = e , leading  to 


y,  A)  dx  dy  dA 


! 


5.4  Simulation 

A standard  Fortran  program  was  written  to  simulate  the  model  described 
in  5.1  and  produce  optimum  density  and  estimates  as  described  in  5.2  and  5.3 
respectively.  The  random  number  generator  used  is  the  same  as  in  the  two- 
dimensional  problem.  This  program  was  debugged  and  executed  on  the  next 
fastest  computer  available  at  NASA  Langley  --CYBER  175,  which  is  faster 
than  the  CDC  6600  by  a factor  of  approximately  3.  Then  a final  program  was  . 
written  for  the  STAR  which  produced  output  which  checked  with  the  C175 
results. 

5.4.1  CYBER  175  Program  --  A point  mass  realization  of  the  conditional 
density  was  simulated  by  choosing  16x96x16  masses  located  inside  a box 
representing  the  three  axes,  the  phase,  phase  rate  and  amplitude,  respectively 
(see  Figure  5.1).  The  16x96  masses  for  the  phase  and  phase  rate  are  located 
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inside  the  rectangle  bounded  by  (-*,*]  In  phase  and  (*  j , £ ] in  phase  rate. 

The  grid  in  the  amplitude  direction  is  floating  according  to  the  location  of  the 
prior  estimate  position.  For  every  new  amplitude  estimate,  the  whole  box  will 
be  centered  on  the  amplitude  estimate  with  eight  masses  equally  distributed 
above  the  center  and  the  same  number  below.  The  mesh  size  between  the 
amplitude  masses  is  chosen  to  be  equal  to  half  the  square  root  of  the 


variance  of  the  estimate. 


Amplitude 
1 K 


^j/A.  


Phase  Rate 
J 


Phase  I ■ 


Figure  5.1 


The  reason  we  chose  the  number  of  discrete  masses  to  be  16x96x16  is  that 
the  maximum  vector  length  for  STAR  operation  is  limited  to  65K  as  well  as 
storage  limitations  on  the  CYBER.  Because  we  know  from  the  two-dimensional 
problem  that  16*96  masses  were  not  accurate  enough  for  density  representation 
the  purpose  of  the  16*96*16  scheme  is  just  to  debug  the  program  on  the  CYBER 
and  try  to  produce  the  same  results  on  the  STAR  later.  This  way,  we  can 
have  a valid  comparison  for  the  speed  of  the  STAR  as  compared  to  other  compu- 
ters for  the  same  amount  of  computation. 

The  main  structure  of  the  program  follows  essentially  the  original 
program  for  the  two-dimensional  phase  demodulator.  The  interpolative  and 
scrambling  matrices  were  precomputed  once  and  for  all,  and  called  out  inside 
the  mass  convolution  loop  for  each  iteration.  Because  of  the  choice  of  the 
dynamic  model  for  the  phase  rate,  a^(y  - z)  is  symmetric  and  constant.  The 
cut-off  limit  for  the  significant  terms  of  a^y-z)  is  determined  by  a radius 


of  length  * Ao  qja  . This  Is  not  the  case  for  the  amplitude  dynamics 
which  does  not  have  a symmetry  because  of  the  damping  factor  and  all  the 
sixteen  terms  have  to  be  computed  for  every  iteration.  When  this  program  was 
executed,  it  produced  estimates  every  12  seconds. 

5.4.2  CPC  STAR-100  Program—  The  standard  FORTRAN  program  in  5.4.1  was  changed 
to  adapt  to  the  vector  notations  for  the  STAR.  The  three-dimensional  array  was 
swept  in  the  right  order  I (phase)  - J (phase  rate)  - K (amplitude).  The 
main  steps  of  the  program  follow  the  same  structure  of  the  two-dimensional 
STAR  program,  except  for  the  added  third  dimension  as  follows: 


(i)  In  order  to  perform  the  scrambling  according  to  the  rule  (I-J)  Inside 
the  convolution  in  the  most  efficient  way,  a block  mapping  VXTOD  is  used  with 
length  17  for  96x16  blocks — see  Figure  5.2.  Of  course,  that  required  an  initial 
step  of  copying  the  box  on  itself  by  using  the'MERGE1  COMMAND  and  carrying  a 
copy  of  the  three  axis  array  for  the  density  masses. 

(11)  The  interpolation  between  two  successive  rows  was  carried  on  the  box  of 
size  17x96x16. 
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(iii)  ‘COMPRESS  COMMAND  was  finally  used  to  get  rid  of  the  extra  17^  ele- 
ment and  to  bring  the  box  into  its  normal  size  16x96x16. 

(iv)  At  this  step  the  box  was  expanded  in  J direction  far  enough  to  cover 
the  maximum  number  of  terms  of  a^y-z)  as  shown  in  Figure  5.3. 


Figure  5.3 

(v)  The  convolution  step  was  executed  inside  a triple  loop  for  the  number 
of  terms  of  a^y-z),  the  third  dimension  k and  the  dummy  variable  for 
integration  over  the  third  dimension.  This  means  that  the  length  of  vector 
for  each  step  was  only  16x96. 

(vi)  Multiplication  by  the  sensor  terms  was  done  in  one  vector  operation  of 
length  24,576. 

(vii)  Normalization  constant  was  obtained  by  folding  the  density  successively 
on  itself  until  it  collapsed  to  single  element. 

(viii)  The  density  then  was  multiplied  by  the  inverse  of  the  constant  and 
transferred  to  the  old  density  for  new  iteration. 


(ix)  This  normalized  density  is  used  for  producing  estimates  for  phase  and 
amplitude.  It  took  89  milliseconds  cpu  time  for  every  iteration  to  produce 
the  same  results  as  on  the  CYBER  175.  This  present  speed  factor  between 
the  STAR  and  CYBER  will  be  cut  down  drastically  after  converting  the 


standard  FORTRAN  program  for  the  CYBER  into  a serial  version  In  an  analogous 
way  to  the  STAR  program  where  all  the  arrays  are  treated  as  single  vectors. 

6.  EXTENSIONS  AND  FUTURE  RESEARCH 

As  is  clear  from  equation  (3.1),  the  computational  problem  of  realization 
of  nonlinear  filters  in  general  consists  of  two  tasks;  first,  the  passage 
from  the  filter  density  to  the  one-step  predictor  density,  a convolution  task; 
and  secondly,  using  the  new  piece  of  data  to  construct  from  the  1-step  predictor 
density  the  next  filter  density.  Independently  of  how  the  densities  are  finitely 
represented,  these  two  partitions  of  the  density  update  are  always  present. 

The  convolution  task  is  time  consuming  when  done  serially  and  calls  for  parallel 
or  pipeline  implementation.  In  general,  strictly  parallel  devices  are 
exhorbitantly  expensive  as  dimension  and/or  griding  increase. 

Currently,  we  are  investigating  optical  realization  of  the  convolution 
task  for  two  state  dimensional  problems  with  Drs.Steier  and  Sawchuck  of  the 
University  of  Southern  California.  The  problem  here  is  one  of  accuracy  and  loop 
closure.  We  hope  to  achieve  with  the  optical  methods  the  accuracy  obtained 
with  the  hybrid  realization.  Another  approach  for  the  convolution  task  would  be 
to  map  space  on  the  line--see  [9]--and  achieve  multi -dimension  convolutions 
by  using  mega-hertz  surface  wave  devices  to  realize  single  dimension  convolutions 
again,  here  the  problem  of  loop  closure  is  important  and  unsolved.  However, 
this  latter  method  is  not  confined  to  two-dimensions  as  the  optical  approach  is. 

Probably  the  most  promising  methods  in  terms  of  high-speed,  low-cost 
systems  are  contemporary  mini -computers  used  with  special  floating  point 
array  processors,  with  the  convolution  task  essentially  done  in  the  array 
processor  and  the  mini -computer  used  to  accomplish  the  second  task. 


An  example  of  such  an  array  processor  with  theoretical  floating  point  arithmetic 
speeds  only  2.78  slower  than  the  STAR  100  is  the  AP-120B  of  Floating  Point 
Systems  Inc.  We  plan  to  program  and  evaluate  such  a system  in  the  course  of 
the  next  year.  Most  probably,  special  purpose  digital  pipeline  devices  offer 
the  most  promise  when  accuracy  is  important,  while  optical  and  other  analog 
devices  would  be  useful  for  filtering  problems  when  accuracy  is  not  that 
important. 

Another  problem  which  merits  attention  is  what  are  the  interactions  of 
other  methods  of  representing  the  density  besides  the  point  mass  used  here,  and 
the  use  of  pipeline  machines.  It  is  not  at  all  clear  that  either  the  Spline  or 
Fourier  filters  can  be  made  to  run  30  times  faster  on  the  STAR  than  on  the  6600; 
essentially  because  the  vectors  involved  are  a factor  of  10  smaller  and  start-up 
times  are  no  longer  negligible. 

Finally,  there  are  two  other  contemporary  machines  which  may  offer 
interesting  possibilities  for  producing  fast  nonlinear  filters;  they  are  the 
Cray  machine  and  a vector  machine  produced  by  Texas  Instruments.  We  hope  to 
be  able  to  report  on  these  in  the  course  of  the  year. 
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APPENDIX 


A.  TWO-DIMENSIONAL  CDC-6600  PROGRAM 


C < 3UCY>S T A F . FOR; 3 4-NOV-76  10:01:44  EDIT  BY  BOCY 

PROGRAM  CYCLIC (INPUT=129,0UTPUT= 129, TAPE5= INPUT, TAPE6*OUTPUT) 

C DESCRIPTION  OF  INPUT  PARAMETERS 

C 

C Y1FST, Y2EST  - THE  EXPECTED  VALUE  OF  INITIAL  POSITION 

C ALP119  - STEADY  STATE  ERROR  VARIANCE  IN  DECIBELS 

C D5LF  - THE  RATIO  OF  DELTA  TO  FILTER  TIME  CONSTANT 

C 222C  - THE  CONTINUOUS  DRIVING  VARIANCE 

C NU  M 1 , NU M2  - ARE  USED  IN  CYCLIC  AND  PROBE  CNLY  AND  COUNT  T 

C NUMBER  OF  PARTITION  POINTS  IN  RECTANGULAR  GRI 

C NO 2 - THE  TOTAL  NUMBER  OF  POINTS  (ESTIMATES)  IN  EACH  SAMP 


DESCRIPTION  OF  DATA  SET 

DATA  MUST  BE  PUNCHED  IN  THE  FOLLOWING  ORDER: 

Y 1 CS", Y TEST, ALP  1 10 , DEL F , Q 22C, NU Ml , NU M 2 , N02 

ALL  REAL  PARAMETERS  (Y1EST  THRU  Q22C)  HAVE  A 10  SPACE  FIE 
ALL  INTEGER  PARAMETERS  (NUN  1 THRU  N02)  HAVE  A 5 SPACE  FIE 
AND  MUST  BE  ?IGHT  JUSTIFIED  IN  THIER  RESPECTIVE  FIELDS. 


CCMME  '•  ?S 


THE  MAI!.’  FLOW  THROUGH  THE  PROGRAM  IS  GOVEE N FD  BY  KOUN'T. 
KOUNT  COUNTS  THE  POINTS  IN  EACH  PATH.  a BLOCK  IS  A SECTI 
OF  THE  PPOGRAM  THAT  HAS  NO  TRANSFER  IN  OR  OUT  EXCEPT 
THROUGH  COMMON. 


COMMON  X DAT  (130,5)  , XHAT(2) 

COMMON  2 (2,  2)  ,PBAR  (2, 2)  ,PN  (2, 2)  ,AN  (2)  ,F  (2,2)  ,PD0HY  (2,2)  , 

♦PDUM Y 2 (2 , 2)  , PNF  (2,2) 

LOGICAL  LOW, UP 

COMMON  /RN/  NZZZ(3),  XNZ(2) 

COMMON  /ON/  DZZZ1,  JGAUSS,  XZZZ  (2) 

COMMON  /PROB/  P 12 , P I , AL P 1 1 0 , DEL F , 22 2C , Y 1 EST , Y 2 EST , 

1 A 11, A2 2, CONST, DELT , FTC, PI DLT , P 1 1 0 , R 1 1 , RX , Q2, Q2  2 
***************************  START  BLOCK  1 ********************** 

2 CONTINUE 
JGA  USS=0 
Y 1 ES7=0 . 0 
Y2EST=0. 0 
ALF1 10=-3.00 
DEI.r=0.  1 
2 22C=  0 . 0 1 
HUM  1*33 
NUM  2=  12  ^ 

N02=1 JO 
N03*  1 

IF (EOF (5) ) 2200,5 
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5 CONTINUE 

WRITE  (6,651)  Y1BST,Y2RST,ALP110,  DBLP,  Q22C,HUH,NHf1?,  *107 
651  FORM  AT (*  *,♦  CYCLIC  INPUT*, 4X,5F10.5, 315) 

P 110  = 10. ** (ALP1 10/1 0.) 

Q 3=0220**  (.25) 

?X=  ( P 1 1 0/  (SQRT  (2.0)  *QQ)  ) **  (4. 0/3.0) 

FTC'^SQRT  (2- 0)  * RX**(.25)  /QQ 
DELT-  DEI.F*FTC 
Q22=  0?2C*PELT 
F 1 1 = RX/D  ELT 

P220=P1 10*S2HT (Q22C/RX) 

IS  AH  P= 1 
NS  AM  F=0 
SUMP  1 = 0. 0 
SUMP  2=0. 0 
CONTINUE 

A 11  = 10.  **( (ALP1 10*1.4) /I  0.) 

A 22=  P220 

K 00  MT= 1 
DELSQ=DELT*+2 
PI=  3.  1415926536 
PI?=2 .0*PI 
PIDLT-=PI/DS  LT 

CONS'r  = -2.')*PiniT’*oTPLr/Q22 
?INV= I.O/^I 
PI2DLT=2. 0*PIDLT 
U 1 = H U M 1 
U2=t::?M2 


XT.I  " = . 75*r*I 
LIMNL=9 
LIME  3 = 0 

v ( 1 , 1 ) = 0 . 0 


r>rc! 

lOi 


3 ( i , 1 ) = o . o 

3 (2, 2)  =022 

A=DELT*PI NV*SQRT (10.0*022) 

I A=  A *0. 5 

IY2  = T’ 2/PI2PLT*SQ?T (50.0*022)  * .5 
CALL  NLF  ( 0, 0 ,Z  1 , Z 2, 5H  AT  , CHAT ,TNLF) 
CALL  PL F (1,0,Z1,Z2, SHAT, CHAT, TNLF) 
X HA T ( 1)  = Y 1 EST 
X H A T f 2)  = Y 2 0 S T 


A=DELT*PINV*SQRT  (10 
I A=  A *0. 5 

IY2  = »,2/ri2PLT*SQ>T( 
CALL  NLF (0,0 ,Z 1 ,Z2, 
CALL  NLF  ( 1 , 0 , Z 1 , Z 2, 
X HA  T ( 1)  = Y 1 EST 
XHAT  (2)  =Y2EST 
X HAT1 =Y 1 EST 
vtiAT2=Y2EST 
eS(1,1)*Al1 
FN  (2,  2)  = A22 
PS  (1, 2)  = 0. 

PN  (2,  1)  =0. 

~ = R 1 1 

d=:vi=  soht(aii) 
DSN=1.0/(?N  (1,1)  *2) 
F (1,  1)=1.0 
F (1,2) =DELT 
F (2, 1) =0. 

F (2, 2)  = 1 . 0 
0EV2=  SQRT  (A22) 


1 


MOD  2 ?I*,2X,*EST.  POSIT.  *f9i. 


DEV 3-  SQRT(RII)  ulji  rnOV' 

CALL  GAUSS(JSEFD,DEV1,Y1EST,X1)  ' »JLL  I 

KOONT=1 

XDAT  (KOUNT,  1)=Xl 

CALL  GAUSS  ( JSEED, DEV 2 , Y2F ST f X2) 

CALL  GAUSS ( JS E ED , DE V 3,  COS { X 1 ) ,Zl) 

CALL  GAUSS  (JSFED,  DEV  3,3  IN  (XI)  ,Z2) 

PEVQ2  r SQ?T(Q2?) 
n=r  1 1 

WRTTF  (6,  1 500) 

FOF'!AT(*0*,3X,*rOSIT.*,5X,*POSIT.  MOD  2 ?I*,2X,*EST.  POSIT. *,9X, 
**71  AMD  7.2*,  1QX,*CYCLTC  LOSS*,5X,*  K-B  EST.  AND  P11*) 
r,n  ’O  4 70 

****************************  END  BLOCK  1 *********************** 

***************************  START  BLOCK  2 ********************** 

CONTINUE 

X1=X1  + X2*DELT 

YDA7 (FOUNT, 1)  = Xl 

CALL  GAUSS(JSEED,DFVQ2, X2,X2) 

CALL  GAUSS  (JS  F ED,  DEV  3 , COS  (XI)  ,Z1) 

CAT!  GAUSS  (JSEED,  DEV  3,  SIN  (XI)  ,Z2) 

***********************  RICCATI  equation  update  **************** 
PDUF.Y  (1  , 1)  = (P*  (PM  (1  , 1)  *2.0*PN  (1, 2) *DELT) -PN  (1 , 2)  **2*DELSQ)  *DEN 
* * PN  (2,2)  ♦DELSQ 

PPUKY  (1,  ?)  =PN  (1,2)  * (R-PN  (1,2)  ♦DELT)  *DEN  ♦ PN(2,2)*DKLT 
PPUBY  (2,2)=-PN  (1,2)  **2*DEN  ♦ PN(2,2)  ♦ Q(2,2) 

PN  (1,  1)  =PDUKY  (1,1) 

PN  (1,2)  = P DH  M Y (1,2) 

PN  (2,2)=  PDTJ  MY  (2,2) 

?N(2,  1)  =PN(1,7) 

DEV  = 1.0/  (PN  (1,1)  ♦ R) 

****************************  i;nD  3LOCF  2 *********************** 

CONTINUE 

'ALL  NLF{1, 1,21,22, SHAT , CHAT , TNLF) 

WRITE  (B,56B7)  TNLF 
FOP*  AT  (FI  0.r>) 

CXHAT  = ATAN2 (SHAT, CHAT) 

PL05S=2. 0*  (1.0-SQRT  (SHA’r**2*CHAT**2) ) 

X DAT ( KOC NT, 2) =CXHAT 
XDAT  (FOUNT,  3)  =PLOSS 
?V?  (1,  1)  = ?N(1,  1)  * 7* DEN 
PNF  (1,2)  = "1N  (1,2)  * 7* PEN 
F N F (2,1)  = PKF  (1,2) 

PNF  (2,2)  =PN  (2,  2)  - PN  (1,2)  **2*DEN 

* ** ****** « ***************** * filter  update  ********************* 

SINF1=SIM (XHA"1) 

COS?  1 = COS  (XHA"’  1) 

X HAT  ( 1)  * X l! A T 1 ♦DEN*  ( - P N (1 , 1)  *STNPl*ZWPN<1,1)  *COSFl*Z2) 

X HAT  (2)  = XHAT20FN*  ( - PH  ( 1 , 2)  * S IN  F 1 *7. 1 *PN  (1 ,2  ) *COSFl*Z2) 

X 1FF=XHA7  (1) 

CONTINUE 

IT (X1FF.GT.PI)  GO  TO  3S 
TF  (FI  FF.  LT. -PI)  GO  T7  U9 

C(3  TO  70 


l “ras-  i'*"  iAwnhwiPA 


:^VV, 


X 1FP=  X 1 FP-PT2 
GO  TO  34 
X IFF*  XI FF*PI2 
GO  TO  «4 
CONTINUE 

IF  ( ARS (CXHAT)  .GT.XLIU) LIN  NL=  LlflNL* 1 
IF  (AUG  (X  1 FF)  .GT.XLIK)  LIN KB=  LI MK B*  1 

*************************  PREDICTOR  UPDATE  ******************* 

X l!  AT  1 * XH  AT  ( 1 ) ♦ PELT  * XH  AT  (2) 

X!1AT2=XHAT (2) 

XDAT  (KOU NT, 4) = XHAT  (1) 

X DAT  ( KOU  NT  , 5)  = PNF  (1,1) 

X 1**OP  = Xl 
CONTINUE 

IF  (XI  MOD.  GT.  PI)  GO  TO  18R 
IF  (XI  NOD. LT. -PI)  GO  TO  1«9 
GO  TO  19  0 
X1N0D=X1  NOD -PI  2 
GO  TO  134 

X1NOD=X1NOD*PI2  T 

GO  TO  194 

CCN’T TIDE  ; '/  , 

ncmifo  n K^/L4/?/  r 

X IF  2-  AES  (THAT  (1 ) - XI ) ‘ / > i;  / 

IF(X1F2.  GT.  PI)  C.l  TO  140  **■  L UPV 

GO  TO  141  V / 

CONTINUE 
X 1F2=X1F2-?I? 

TKFUI.r  = 1K35LP*1 
go  to  in 
CONTINUE 

E?»>LF  = ABS  (XI  FF-X1NOD) 

ERF NL= AOS  (CX MAT- XI TOO) 

IF  (SPPLF.  GT.  P I ) LP  PL  F = A D G (EPRIF  -PI2) 

IFfErPVL.GT.  FI)  ERRf'L*  AI3G  (FRFNL-PI2) 

SFRJF* ABS  (X1MOD-CXHAI) 

K ? IT  (R  ,201)  FOUNT,  XDAT  (FOUNT,  1)  , X 1 KOD  , XDA  T ( * CUNT  , 2)  ,21  ,22,  (XDAT 
* (KCtlNT,!)  ,1=3,5) 

F0  7.1  (*0*,  13,  IX,  1P3E14.A,4X,  1?2E14.^,4X,  1P1E  14.R  /) 

IF  (EOnjJT.TQ.  N02)  GO  TO  505 

KOU  NT  = KOU  NT  ♦ 1 

GO  TO  410 

CONTINUE 

GUFP=0.0 

S U fl  C = 0 . P 

DO  1501  1=31 , NO 2 

X D=  ABS (XDAT (I, 1)  -XDAT  (1,2) ) 

CONTINUE 

IF(XD.GT. PI)  GO  TC  1499 

GO  "*0  1500 

XD=XO-PI2 

GO  TO  1498 

StJNP=  (XU)  **2  ♦ SUXP 

S UKC*  X DAT  (I, 3) ♦SUEC 

CONTINUE 


mew* 


■(PjMqwWirHJUffn  •' 


r 


r 


H=N02-10 

SUHP=SIINP/H 

sunc= sumc/h 

XNSAMP^NSANP 

XAA=XNSANP*1.0 


SUMP1*  (SUMP*XNSANP*SUNP1)  /X  A A 
dsunm  = alooio  (sumpI)  *io. 

KPTTF  (6,  150  3)  • 

1508  FCRNA','(*0*,5X,*NONLINSAr  CYCLIC  ESTIMATOR*) 

WRITE  (6  , 1511)SUHP1,  DSUMP1 

1511  FOPMA?  (*0*,*  AV  ER  AGE  STATISTICAL  VARIANCE  =*,1PE13.6,  10X, 

* * AV  TPAGE  COMPUTED  VARIANCE  = *, 1 PEI  3. 6//) 

SUMP=0. 0 
S»MC=O.0 

DO  1601  1=31, NO? 

XD= ABS  (XDAT  (1,1) -XL AT  (1,4) ) 

16  93  CONTINUE 

IF(XD.GT.PI)  GO  "O  16  99 
GO  TO  1700 
699  XD=XD-?I2 
GO  TO  1698 

1700  SUMP=  (XD)  **2*SUMD 

S U MC=  XDAT  (1,5)  ♦StiMC 
If- Cl  CONTINUE 
S DM  P=  SUM  P/H 
Side's  SUMC/H 

5UMP2  = (SUMP*XNSAMP*S0M?2) /XA A 
DSUMP2=ALOG10 (SUMP?) *10. 

WRITE  (6, 1509) 

SOQ  FORMAT(*9*,5X,  *?E-5i:JEA  HIZED  K-3  FIL'E?*) 

WRITE  (6,1511)  SUFP2,  D3UMP2 
NSAf!P=!iS  AMP*  1 

TF(ISAM'>.  S3.  NO 3)  GO  TO  2200 
ISA  MP  = I SAMP* 1 

Go  to  11 

****************************  BLOCK  1 *********************** 

1 20  0 VfRT?P  (6,2201) 

’201  FCP/A'r  (*D*,99j:,*vopMAL  COMPLETION*) 

STOP 

END 


$*. 

\ $ 


& 
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r 1 

3UBFOU7INE  GAUSS(JS,SD,XH,X) 
DIMENSION  NST (2) 

COMMON  /RN/  HI,  N2,  HC,  T1,  T2 
COMMON  /GN/  TWOPI,  J,  XR  (2) 

TF  (U)  10,  10,  20 

10 

J = 2 

oy.\ 

* 

TWOPI =3. *ATAN (1 .) 

NST  (1) =102943 

NST (?) = 1955 17 

• 

XP  ( 1) =3ANF (NST, 1) 

GO  ”0  35 

/y// 

«>  f Jr  A 

20 

GO  TO  (30,40) , J 

30 

J = 2 

XP  (1)  = 13 AN F (NST , 0) 

% 

35 

X R ( 2)  =DANF  (NST,  0) 

X 1=SgRT  (ARS  (-2.  *ALOG  (XR  (1)  ) ) ) 
XR  (2)  =TWOPI*XR  (2) 

XP  (1)  =X1  *SIN  (XR  (2)  ) 

XR  (2)  =X1*C0S  (XP  (2)  ) 

X = XR  (1) *SD»XM 

RETURN 

4 0 

J = 1 

X = XP  (2)  *SDO:M 

IF?  UP  N 

FND 

r» 


10 


100 

no 

120 

100 

200 


T1.  T2 


FUNCTION  3A  VF  (NS,  MODE) 

DIMENSION  US  (2),  NC  (2) 

COMMON  /F<  M/  Nl,  U2,  UP  , i., 

DATA  Ml,  M2/244'»34,  15355  1/ 

iode=o  ra  continue,  otherwise  pestapt 

INTEGER  NUMBER  NS  ( 1)  *2**  13*NS  (2) 

IF  (MODE)  10,  100,  10 
N1=NS  (1) 

N2=N’S  (2) 

T 1=  ?.  **  (-  1R) 

T 2=2 . ** (-30) 
l«)P=2**1  o 
00  200  1=1,2 
GO  TO  (1  10,120)  ,1 
K=M2*N2 
GO  TO  190 
K=Hl*N2*M2*NH!<P 
KD= K/MP 
NC  (I)  =F-KD*rr 
N 1=  NC  (2) 

N2=»'C  (1) 

X N 1 = N 1 
v n2=  V2 

EAN7=XNl*T1^XS2*T? 

PETES  N 
END 


WITH 
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SUBROUTINE  NLF  (HC,S AHP, Z 1 , Z2 , S HAT r CHAT, TNLF) 

INTEGER  NC,SAHP 

REAL  Z1,Z2, SHAT, CHAT, TNLF 

INTEGER  1,11, JC,J1, J2,K1 ,K2,KL,KH, NJ1 , NJ2 , N K 1 , NK2 , NTERH, 

1 MTERH32,NSI7,E,NSIZE32 

PEAL  All,  A22,Ct,CNOPH,  CONST,  CP,  PI,PIDET.  ,Q22  ,SI,T,TT, 

1 YlFST,Y2EST,TEHP,Tfc!1P1,TENP2 
PfcAT.  IN  (8) 

PEAL  TROtf  (32) 

REAL  COS  Y ( 32)  , SINY  (32)  , SN1  (12)  ,S1  (32)  ,S2  (32 ) , SIGH  A ( 32) 
RF.AL  P3I  (128)  ,A  (10)  , DELJ  (123) 

INTEGFD  JNS  (4224) 

PEAL  JN  (4046)  , JK1  (4096)  , JO  (4  096)  ,JNA  (4  756) 

COHNON  /P ROB/  THOPI, PI, ALP1 10 , PELF, Q22C, Yl EST ,Y2EST , 

1 All, A 2 2, CONST, DEL, FTC, PI  DEL, Pi  1 0 ,RDEL,  R Y , QQ , Q2 2 

COM  HON  /NLFC/  NC, NT , NTERH , NTERH3 3 ,S 1 , S2 ,S IG HA , PSI , A, COS Y , 
1 DELJ,JO,JNA, JNS, SINY 
EQOIVALENCF  (JN  1 (1 ) , JNA  (32  1 ) ) 

IF  (SAMP.LE.O)  GO  TO  100 
SET  CLOCK 
T=  S ECON D (T) 

EVA  LUA^E  5ENSOF  7EFM3 

ro  io  t=i,32 

ir  sni  (i)  = sxp(zi*s  i (i)  ♦z2*sr  (i) ) 

FC'’*  THE  INTERPOLATED  .T  N AND  T'UT  IN  JN  1 
J1  = 0 
J2  = 0 

DO  30  T=  1 , 1 2 8 
TEMP=DELJ  (I) 

DO  20  .1=1,32 
K 1= J 1 ♦ J 
KL= JNS  (K  1) 

K H = J N S (KlO) 

TFMP1=JN  (KL) 

K2=J2*T 

20  JN  1 (K2)  = TSHP  1 ♦TEMP*  (JN(KH)  -TEMPI) 

J1  = JU33 
1C  J 2= J2 ♦ 32 

EXPAND  INTEPPOLATED  HATSTY  ON  BOTH  SIDES 
J 1=  NJ  1 
J2=NJ2 
v 1 = NK  1 
K 2 = V K 2 

DO  40  1=1, MTFPM  32 
JNA  (J  1)  = JNA  (J2) 

J N A (FT  1)  - TNA  (K2) 

J1  = JU1 
J 2= J2* 1 


K 1 r K 1 ♦ 1 


4 0 Y 2=K2* 1 


r, 

i I 

t 


i 

< 

* 


J1»WSISE32 
J2= J 1 

DO  f>0  I = 1 , NTERM 
J 1*J  1*32 
J2=.12-l:> 

TEMP=  A (I) 

DO  60  J=  1,4  096 
K 1=J1*J 
K2=J2»J 

60  JN  (J)  - J N (J)  ♦TEMP*  (JNA(Kl)  »JNA  (K  2)  ) 

C CUMULATE  ROW  SUMS 

no  99  T = 1,32  - _ J 1 i / 

I 1=1 

TEMP2  =JN  ( T 1 ) 

DO  TO  0=1,127 
1 1=  T 1 ♦ 32 

70  TEM?2=TEMP2*JN(T1) 

PC-  TPOW  (T)  = TEM”2 

C ACCUMULATE  ESTIMATES  AND  NORMALIZATION  CONSTANT 

CNOSS=TFOW(1)  *S  N 1 (1) 

SHAT=STNY  (1)  *CNOPM 
CHAT=COSY  (1)  *CNO»r 
DO  95  1=2,32 
TIM  ?2='T,ROW  (T)  *SN1  (I) 

S H A T=  SH  AT *S  I N Y ( I ) *TEHP2 
OHAT  = CHATfCOSY  (T)  *TFI1P? 

85  CNORW=CNORmTEMP2 

C 1 - O/CNOpr 

c yr  7=  SHAT*CNO~.Y 
C:;.AT=CHAT*CN’ORM 

C TRANSFER  NORMALIZED  DENSITY 

DO  ap  7=1,32 

n = r 

Trr.P2  = s;n  (i)  *o?to?.r 

no  =»c  J= 1 , 129 

J‘.  (T  1 ) = TE!SP2*  J V (T  1 ) 

O'*  11  = 11*32 
C 'rT*'"C'0T 

T SI  F = SECOND  (TT)  -T 
RETURN 

c i i " r a l t z e sample  path  by  transfertng  jo  to  jn 

100  IF  (XC.LE.O)  GO  TO  200 
DC  110  1=1,4096 
11C  JN  (I)  =.70  (I) 

T ET  U P V 

C GI.ORAI  INITIALIZATIONS  FOR  NONLINEAR  riLTEF 

200  NSIZE=10 

NTFPM=64 . 0*SQRT (50. *Q22) /PTDEI.*0.5 

IF  (NTEFM.GT.  NSIZE)  N','ERM=NSIZE 

NSTZF32=NSIZE*32 

NTEP*32=NTE?M*32 

NK2=NSIZE32*1 

NJ1=NK2-NTE?N32 

NJ2=NSIZE32*4097-NTEPM32 

NY1  = VSI7, 232  *4007 


t 
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PHASE  VARIABLES 
DO  210  T=1,3? 

SIGMA  (I)=PI*((2.*I-1.)/32.-1.) 

COSY  (I)  =COS  (SIGMA  (I)  ) 

SINY  (T)  =SIN  (SIGMA  (I)  ) 

SI  (T)  -COSY  (I)  /FDEL 
2 10  S2  (T)  = 5 T NY  (I)  /RDEL 

PHASE  HATE  VARIABLES 
DO  220  1=1,123 

220  PST  (T) =PIDEL* ( (2. *1-1.) /I 2B. -1.) 
SETUP  THE  TRANSFER  MATRIX 
DO  240  0=1,123 
J1  = (.1-1)  *32 

0 2r  (0-1)  *33 
DO  230  1=1,32 

11  = 01 *1* MOD (4  6-  (J- 1 ) /4* 1,32) 
12=02*1 

230  0 NS  (12)  =11 

240  JNS  (02*33)=JNS (02*1) 

SETUP  THE  INTERPOLATION  VECTOR 
IN  ( 1)  =0  . 375 
T N ( 2 ) =0.623 
IN  (3)  =0. 3 75 
IN  (4) =0.  1 25 
in  (5)  =in  (i) 

TN(f)=IN(2)  Ki-Y 

1 N ( 7 ) =T  N ( 3 ) yw 

IN  (3)  =TN  (4) 

0 = MOD  (STERN  ,4) 

DO  245  1=1,4 
243  D Z L 0 (T)  =TN  (1) 

DO  :c0  T=5, 125,4 
'OSLO  (T)  =DEL0  (I  -4) 

OSLO  (1*1) =DELJ  (1-3) 

DEI  0(1*2)  =DEL.T  (1-2) 

250  2 EL  J (1*3)  =D  E L J (1-1) 

EVATUATE  CONVOLUTION  TERMS  A (T ) 

PD  2 30  1 = 1 ,NTETM 
TEMr,=  T/1  23. 

TFw?=CONST*TEMP*TF.M  P 
A ( I ) s ? . 

IF  rTr?.GT.-47)  A (T)  =2XP  (TEMP) 

2.30  CONTINUE 

CISSTSOCT  THE  A PRIORI  DENSITY 
CN05r  = 1 . 0/(TMOPI*SQ3T (A  1 1*A22) ) 

C L = - ^ . 5/A22 
~T  = - ' . 5/A  1 1 
DO  230  1=1, 32 
C?='IG*A  (I) -Y1 EST 
CF=C=*C?*SI 
01=0 

DO  230  0=1, 123 
J 2=  H*I 

TEMP=PSI  (0)  -Y2EST 

JO  : T ; ) = EX  P (TEMP*TEM  P*CL  *CP)  *CNOPM 
230  01=01*32 
r rr  <’  r v 
END 
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B.  TWO-DIMENSIONAL  CDC-STAR-lOO  PROGRAM 


’i 


* 

'if 

i 


PROGRAM  MAIN( I NPUT , OUT PUT , T APE5 = I NPUT , T A PE 6=0UT PUT ) 

C SAMPLE  PATH  VARIABLES 

REAL  SXI ( 130 ) » SCHAT ( 130) ,SSHAT( 1 30 ) , SX 1HATNL ( 130 ) , SERRNL (130), 

1 TNLEI 130) » SPLOSSNL (130) ,SX1HATPL ( 1 30 ) , SERRPL ( 1 30 > , S P 1 1 PL ( 1 30 ) 

C f OMUL  A T I VE  SAMPLE  PATH  VARIABLES 

REAL  CERRNU  130)  ,CESQNL(  130  ) ,CEVARNL  ( 130 ) ,C0BNL( 130)  ,CCSNL(  130) 
REAL  CERRPL (130) ,CESQPL( 130) ,CEVARPL( 130  ) ,CDBPL ( 130 ) ,CCSPl ( 130 ) 
C MONTE  CARLO  SUMMARY  STATISTICS 

REAL  XERRNL,XESQNL,XEVARNL,XDBNL,XCSNL, 

1 XERRPL,XESQPL,XEVARPL, XOBPL,XCSPL 

C SINGLE  SAMPLE  VARIABLES 

REAL  CHAT,SHAT,X1HAT,P11,X1  ,Z1  ,Z2 

C constants 

REAL  CS( 130) ,DBEPS,EPS( 130) ,TWOP I , ZERO ( 1 30 ) ,0NE ( 1 30 ) , P 1 2 ( 130) 

C WORKING  VARIABLES 

INTEGER  I,J,K,L 
REAL  T * TEMP ( 130 ) 

LOGICAL  PATH ,CUMPATH 
BIT  BT ( 130) 

C PROBLEM  SETUP  VARIABLES 

REAL  ALP110,DELF,Q22C,Y1EST.Y2EST 
INTEGER  NMC,NSAMP,MD,NO 
C DERIVED  PROBLEM  CONSTANTS 

REAL  All ,A22, CONST tDELt FTC t PI ♦PIDEL,P110,RDEL,RX,00t022 
C PROBLEM  COMMON 

COMMON  /PROB/  TWOPI ,PI , AL  P 1 U)  ,DE LF  , 022C  , Y 1 E ST  , Y2E  ST  , A 1 1 , A 22  , 

1 CONST, DEL, FTC , PI  DEL, PI 10,RDEL,RX, 00,022, MO ,ND 

C SET  PRINTOUT  CONTROL 

PATH=. TRUE. 

CUMPATH=.TRUE. 

C READ  INPUT  PARAMETERS 

10  READ  ( 5,5000,EN0  = 500)  Y1 E ST , Y 2EST , AL P 1 10 , OE L F , 02 2C  ,NMC , NS AMP  ,MD 
5000  FORMAT ( 5E10.4,3I 5) 

C COMPUTE  THE  CONSTANTS 

MD=(MD/2 )*2 
IF  (MD.LT.20)  M0  = 3 2 
IF  (MD.GT.64)  MD=3 2 

nd=<+*md 

PI=4.0*ATAN( 1.0) 

TWOP I =2 .0*P I 
P110=10.**( ALP110/10. ) 

QQ=Q22C** ( .25) 

RX= ( PI  10/ ( SORT ( 2.0 ) *Q0 ) ) ** ( 4 .0/3 .0 ) 

FTC=SQRT (2.0)*RX**( . 251/00 

DEL=DELF*ETC 

Q22=Q22C*DEL 

P I DEL=P I /DEL 

RDEL=RX/DEL 

A1 1 = 10.0** ( ( ALP1 10+1 .4 ) /10.  ) 

A22=2.0*P110/(FTC*FTC) 

CONST=-2.0*PIDEL*PIDEL/022 
CS( 1 : 130 )=0.75*PI 
P I 2 ( 1;130)=TW0PI 
Z ERO ( 1,130) =0. 

EPS( l;130)=l.E-50 
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D6EPS=AL0G10(EPS( 1 > ) 

ONE ( 1 ; 1 30 ) * 1 . 

C INITIALIZE  CUMULATIVE  SAMPLE  PATH  VARIABLES 

CERRNLI 1 ? 130)=0. 

CFSONLI I 5 1 30  ) =0 . 

LCSNL ( 1 5 130) =0. 

CERRPLI 1 ;130)=0. 

CESQPL ( 1 5 1 30 ) =0 . 

CCSPL ( 1 ; 130 )=0. 

c PRINTOUT  PROBLEM  PARAMETERS 

WRITE  (6,6000)  NMC , N S AM P , ALP  1 1 0 , DE LF , 02 2C , Y 1 E ST , Y2E ST , 

1  P110,QQ,RX,FTC»DEL»Q22»PIDEL,RDFL»A11,A22»C0NST»CS( 1 ) 

6000  FORMAT ( 1H1,31X, 18HPR0BLEM  PARAMETERS/ 

1 1H0,11X,9HPARAMETER,6X,5HVALUE , 1 1 X ,9HPARAMFTER  , 

2 6X, 5H VALUE/ 1H0, 14X , 3HNMC , 9 X , 14 , 1 4X , 5HNS AM P , 8 X , 14/ 

3 13X,6HALP110,3X,E15.8,8X,4HDELF,4X,E15.8/ 

4 14X,4H022C*4X,E15*8*8X,5HY1EST,3X,E15.8/ 

5 14X,5HY2EST,3X,E15.8,8X,4HP110,4X,F15.8/ 

6 15X,2HQ0,5X,E15.8,9X,2HRX,5X,E15.8/ 

7 15X,3HFTC,4X,E15.8,9X,3H0EL,4X,E16.8/ 

8 15X,3H022,4X,E  15. 8 , 8X , 5HP I OE L , 3X , F 1 S . 8/ 

9 14X,4HR0EL,4X,E15.8,9X,3HA11 ,4X,E15.8/ 

A 15X,3HA22,4X,E15.8,8X,5HCONST,3X,Fl‘>.8/ 

B 15X,2HCS,5X,E15.8) 

C BEGIN  THE  MONTE  CARLO  PROCESS 

C INITIALIZATION  OF  THE  SUBROUTINES 

CALL  STATE ( 0,0, XI ,Z1 ,Z2 ) 

CALL  NLF(0, 0,0.,0., SHAT, CHAT, T) 

C MONTE  CARLO  LOOP 

DO  200  K = 1 , NMC 

C INITIALIZATION  OF  SAMPLE  PATH  VARIABLES 

CALL  STATE(K,0,X1,Z1,Z2) 

CALL  NLF(K,0,Z1,Z2, SHAT, CHAT, T) 

CALL  PLL(K,0,Z1,Z2,X1HAT,P11 ) 

C SAMPLE  PATH  LOOP 

DO  100  J=1 , NS AMP 
CALL  STATE(K,J,X1,Z1,Z2) 

CALL  NLF(K, J,Z1,Z2, SHAT, CHAT, T) 

CALL  PLL(K,J,Z1,Z2,X1HAT,P11  ) 

C STORF  THE  SAMPLE  VARIABLES 

SX1 ( J )=X1 
SSHAT( J)=SHAT 
SCHATt J)=CHAT 
TNLF ( J ) =T 
SX 1HATPL ( J ) =X 1HAT 
100  SP11PLI J ) =P  1 1 

C VECTOR  ACCUMULATE  THE  S AM  PL  F PATH  AVERAGES 

SX1HATNL ( 1 ; 130)=VATAN2( SSHAT ( 1 ; 130) ,SCHAT( 1:130); 

1 SX 1 H ATNL ( 1 1 1 30 ) ) 

SERRNL ( 1 5 130 ) =SX 1 ( 1 : 1 30  I -SX1H ATNL ( 1 J 1 30  ) 

CALL  M0D2PI ( SERRNL) 

SPLOSSNL ( 1 5 130) =SSHAT (1;130>*SSHAT(1:130)+SCHAT(1:130)* 

1 SCHAT (1*130) 

SPLOSSNL ( 1 ; 130) =2.0* ( 1 . O-VSOR T ( SPLOSSNL (1:130); SPLOSSNL (1:130))) 
CERRNL( 1 5 130 ) =VAVG(CERRNL (1:1 30 ) , SERRNL ( 1:130) ,K :CERRNL (1:130)) 
.TEMPI  1 ; 130)=SFRRNL( 1 : 130)*SERRNL (1:130) 


' FPn 


— — ‘ 


CESONL  ( 1;  130)  =VAVG(CESQNL  ( 1 ; 130 ) , TEMP  ( 1 ; 130  ) , K ; CE  SQNL  ( 1 ; 1 30  ) > 

IF  (K.LE.l)  GO  TO  110 

CEVARNL < 1 ; 130  )=CFSONL  ( 1 ; 130 )-CFRRNL  < 1 ; 130)*CERRNL  (1:130) 

H T < 1 ; 1 30 ) = ( CEVARNL (1:130). LE.EPS(1;130)> 

CEVARNL  (1:130)  =QRVMASK( EPS(  1:130)  t CEVARNL  ( 1 ; 1 30  ) , BT  ( 1 : 1 30  ) : 

I CEVARNL ( 1 : 130)  ) 

l.ORNL (1:130) =VAL0G10( CEVARNL (1:130) ;CORNL (1:130)) 

CDRNL ( 1 : 130)=10.0*C0HNL (1:130) 

110  HT ( 1 ; 130) =(SERRNL( 1 ; 1 30 ) ,GT .CS ( 1 5 1 30 ) ) 

TEMP ( 1 ; 130)=Q8VMASK( ONE (1:130), ZFRO ( 1 ; 1 30 ) , BT (1:130): 

1 TEMP(l;130)) 

CCSNL (1:1 30)=VAVG(CCSNL(  1;130),TFMP(1;130) ,K:CCSNL (1:130)) 
SERRPL ( 1 ; 1 30 ) =SX 1 ( 1 ; 130 ) -SX 1 HAT  PL (1:130) 

CALL  M0D2PI ( SERRPL ) 

CFRRPL ( 1 : 130 )=VAVG(CERRPL ( 1 : 130 ) , SERRPL ( 1 ; 1 30 ) , K ; CFRR PL ( 1 ; 1 30 ) ) 
TEMP ( 1 ; 1 30 ) =SERRPL ( 1 ; 130)*SFRRPL( 1:130) 

CESOPL ( 1 ; 130)=VAVG(CES0PL ( 1 : 130) , TEMP ( 1 : 1 30 ) , K ; CE SOPL (1:130)) 

IF  (K.LE.l)  GO  TO  120 

CEVARPLf 1;130)=CESQPL( 1:130)-CFRRPL( 1:130) *CEPR PL (1:130) 

4T(  1 : 130)  = (CEVARPI.  ( 1 ; 130)  ,LF.FPS(  1;  130)  ) 

CEVARPL ( 1 ; 130) =08 VM ASK (EPS ( 1 : 1 3 0 ) , CE V A« PL ( 1 : 1 30 ) , HT ( 1 : 1 3 0 ) : 

1 CEVARPL ( 1 ; 130) ) 

COBPL ( 1 : 130 )=VAL0G10( CEVARPL ( 1 : 1 30 ) : COR PL ( 1 : 1 30 ) ) 

CDBPL ( 1 : 130)=10.0*COBPL( 1 J 130) 

120  HT ( 1 ; 1 30 ) = ( S ERR P L ( 1: 130).GT.CS( 1:130)  » 

TEMP ( 1 ; 130 ) =08 VM ASK ( ON 6 ( 1 5 1 30 ) , ZFRO ( 1 : 1 30 ) , HT ( 1 : 1 30 ) : 

1 TEMPI  1:130) ) 

CCSPLI 1 ; 130)=VAVG(CCSPL ( l: 130) ,TFMP( 1 : 130  ) ,K  :CCSPL  ( 1 : 130)  ) 

IF  ( PATH. AMD. (K.FO.l  ))  GO  TP  130 
GO  TO  200 

C PRINT  SAMPLE  PATH  VARIABLES 

130  WRITE  (6,6010 

6010  FORMAT! 1H 1 , 36X , 2 1HS AMPLE  PATH  VARIABLES/ 

1 38X  , 1 9H ( FIRST  SAMPLE  PATH)/) 

WRITE  (6,6011) 

6011  FORMAT ( /18X.9HSX1-  ,6X, 

1 36H  PHASF  VARIABLES  / 

2 20X  , 1HI  , 13X  . 1 1HSX1 ( I ) ,8X, 13HSX1 ( 1 + 1 ) /) 

WRI TF< 6,6100) ( ( I , SX1 ( I ) ,SX1( 1+1 ) ) , 1=1, 129,2 ) 

WRITF  (6,6012) 

6012  FORMAT ( /18X,9HSSHAT“  ,6X, 

1 36H  SIN(SXl)  ESTIMATES  / 

? 20X , 1HI , 13X, 1 1HSSHAT ( I ) , 8 X , 1 3H S SH AT ( I + 1 > /) 

WRI TE(6,61C0)  ( ( I ,SSHAT ( I ) ,SSHAT ( I + 1 ) ) , I = 1 , 1 29 , 2 ) 

WRITE  (6,6013) 

6013  FORMAT  ( /18X,9HSCHA'f-  ,6X, 

1 36H  COS(SXl)  FSTIMATES  / 

2 20X,1HI , 13X , 1 1HSCHAT ( I ) , 8X  , 1 3H SCH AT ( I + 1 ) /) 

WRI TE (6,6100) ( ( I, SCH AT ( I ),SCHAT(  I + 1 ) ) , I = 1 , 1 29 , 2 ) 

WRITF  (6,6014) 

601^  FORMAT ( /18X,9HSX 1HATNL-  ,6X, 

1 36H  NONLINEAR  ESTIMATES  / 

2 20X,1HI ,13X, 11HSX1HATNL( I ) , RX , 1 3HS X 1 H ATNL ( I +1 ) /) 

WRITE (6,6 100)  ( ( I »SX  1HATNL ( I ) ,SX1HATNL<  I + 1)),I  = 1,129,2) 

WRITE  (6,6015) 

6015  FORMAT! /18X,9HSERRNL-  ,6X, 

l 36H  NONLINEAR  ERRORS  / 
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J 1 J J '**J  L.L. 


2 20X , 1HI ,13X,11HSERRNL( I ) , 8 X , 1 3H  SER  RNL  ( I + 1 ) 

WRI  TE (6,6100)  < ( I ,SERRNL ( I ) , SFRRNL (1+1)1,1=1,129,2) 
WRITE  (6,6016) 

6016  FORMAT ( /18X,9HSPL0SSMl-  ,6X, 

1 36H  NONLINEAR  RLOSS  FUNCTION  / 

2 20X,1HI ,13X,1 1HSPL0SSNH I ) , 8 X , 1 3HSPL0SSNL ( 1 + 1 ) 

WRI TF( 6,6100) (< I , SPLOSSNL ( I ) , SPLOSSNL ( I +1 ) ) , I = 1 , 1 29 , 2 ) 
WRITE  (6,6017) 

6017  FORMAT ( / 1 8 X ,9HTNLF-  ,6X, 

1 36H  NONLINEAR  EXECUTION  TIMES  / 

2 20X , 1HI , 13X,  1 1HTNLF ( I ) , 8X , 1 3HTNL F ( I + 1 ) /) 

WRITE(6,6100)  ( ( I ,TNLF(  I ) ,TNLF(  1 + 1 ) ) , 1 = 1,129,2) 

WRITE  (6,6018) 

6018  FORMAT ( /18X,9HSX1HATPL—  ,6X, 

1 3 6 EH  PHASE  LOCK  LOOP  ESTIMATES  / 

2 20X,1HI ,13X,11HSX1HATPL(  1 ' , 8X , 1 3H S X 1 H AT  PL ( I + 1 ) 

WRITE(6,6100)((I,SX1HATPL(I),SX1HATPL(I+1)),I=1,129,2) 
WRITE  (6,6019) 

6019  FORMAT) /18X,9HSFRRPL-  ,6X, 

1 36H  PHASE  LOCK  LOOP  ERRORS  / 

2 20X.  1HI  ,13X,  llHSFRRPLf I ) ,8X, 13HSERRPL ( 1 + 1 ) 

WRITE (6, 6 100) ( ( I ,SERRPL( I ) , SERRPLI 1 + 1 ) ) , 1 = 1 ,129,2  ) 
WRITE  ( 6,6020) 

6020  FORMAT ( /18X,9HSP1 1 PL  — ,6X, 

1 36HP11  FROM  PHASE  LOCK  RICATTI  EQUATION  / 

2 20X,1HI , 13X  , 1 1HSP1  1PL ( I ) ,«X,13HSP11PL(  T + l ) 

WRI  TF(6-,6100)  ( ( I , SP1  13L  ( I ) ,SP1  1PL(  1 + 1 ) ) , 1 = 1 , 120,2  ) 

200  CONTINUE 

PRINT  CUMULATIVE  SAMPLE  PATh  VART  A8LP.S 
IF  ( CUMPATH.AM'-).  ( NMC.GT  . 1 ) ) OP  TO  310 

c-n  to  400 
310  WRITE  (6,6030) 

6030  FORMAT(  1H1,30X,32HCUM()LATIVE  SAMPLE  PATH  VARIAKLFS  /) 
WRITE  (6,6031) 

6031  FORMAT ( /18X ,9HCERRNL-  ,6X, 

1 36W  CUMt  <L  AT  I VP  NONLINEAR  ERRORS  / 

2 POX , 1HI  , 13X , 1 1HCESRNL ( I ) , 8X , 1 3HCFR RNL ( I + 1 ) /> 

WRITE  (6,6100  )((  I , C F R R.  N L ( I ) * C F R R n L ( I + 1 ) ) , I = 1 , 1 29 , 2 ) 
WRITE  (6,6032) 

6032  FORMAT (/ 18X ,9HCERRPL-  ,6X, 

1 36HCHMULATI VE  PHASE  LOCK  ERRORS  / 

2 20X,  1HI ,13X,11HCERRPL(  I ) , 8X , 1 3HCE RR PL ( I + 1 ) /) 

WRITE  (6,6100 ) ( ( I ,CERRPL ( I ) ,CERRPL ( 1 + 1 ) ) , 1 = 1 ,1 29,2  ) 
WRITF  (6,6033) 

6033  FORMAT ( /18X,9HCES0NL-  ,6X, 

1.  36HCUMUL  AT  I VE  NONLINEAR  SOU  ARE  0 ERRURS  / 

2 20X,1HI ,13X,11HCFS0NL( I ) ,8X , 13HCES0NL ( I +1 ) /) 

WRITF  (6,6100) (( I , CF  SONL ( I ) ,CFSONL ( I + 1 ) ) , I = 1 , 1 29 , 2 ) 
WRITF  (6, 603 A) 

6034  FORMAT ( /18X,9HCESQPL-  ,6X, 

1 36HCUMULAT I VE  PHASF  LnCK  SQUARED  ERRORS  / 

2 20X, 1HI , 13X, 1 1HCESQPL ( I ) , 8X , 1 3HCF SQPL ( I + 1 ) /) 

WR ITF  (6,6100  )((  I , C F S 0 P L ( I ) , CF SQPL  < I + 1 ) ) , I = 1 , 1 29 , 2 ) 
WRITF  (6,6035) 


6035  FORMAT (/18X,9HCE' - RNL-  ,6X, 

1 36HCUMULATI VE  NONLINEAR  ERROR  VARIANCE  / 

2 20X,1HI ,13X,11HCEVARNL(  I ) , 8X , 1 3HC E V ARNL II  + 1 ) /) 

WRITE  (6,6100) ( ( I .CE  V ARNL ( I ) , CE VARNL ( 1 + 1)1,1  = 1,129,2) 
WRITE  (6,6036) 

6036  FORMAT ( / 18X ,9HCEVARPL-  ,6X, 

1 36HCOMULAT 1 VE  PH ASF  LOCK  FRROR  VAR1ANCF  / 

2 20X , 1HI , 13X, 1 1HCEVARPL ( I ) , 8X , 1 3HCE V AR PL ( I + 1 ) /) 

WRITE  (6,6100)  ( ( I ,CFVARPl ( I ) ,CEVARPl ( 1 + 1 ) ) , 1 = 1 , 120,2 ) 
WRITE  (6,6037) 

6037  FORMAT ( /18X .9HCDBNL-  ,6X, 

l 36HCUMUL AT  I VE  NONLINEAR  ERROR  DECIBELS  / 

? ?0X,1HI ,13X,11HCDBNL( I ) , RX , 1 3HCDBNL ( I + 1 ) /) 

WRITE  (6,6100) ( ( I ,CDBNL( I ) ,COBNL( 1 + 1 ) ) , 1=1 ,129,2  ) 
WRITE  (6,6038) 

6038  FORMAT ( /18X,9HCDBPL-  ,6X, 

1 36HC  OMUL AT  I VE  PHASE  LOCK  ERROR  DECIBELS  / 

2 20X, 1HI , 13X, 1 1HCDBPL( I ) , 8X  , 1 3HCD8PL ( I + 1 ) /> 

WRITE  (6,6100  ) ( ( I ,CDRPL ( I ) ,COBPL(  1 + 1 ) ) , 1 = 1 ,129,2 ) 
WRITE  (6,6039) 

6039  FORMAT ( /18X,9HCCSNL~  ,6X, 

1 36HCDMULATI VE  NONLINEAR  CYCLE  SLIPS  / 

2 20X,1HI ,13X,11HCCSNL( I ) , RX  , 1 3HCC SNL ( I + 1 ) / /) 

WRITE  (6,6100) (( I ,CCSNL( I ) ,CCSNL ( 1+1 )), 1=1 , 129,2 ) 
WRITE  (6,6040) 

40-^0  -ORMAT  ( / 18X  ,9HCCSPL-  ,6X, 

1 36HCDMUL AT  I VE  PHASE  LOCK  CYCLE  SLIPS  / 

? 20X , 1HI , 13X, 1 1HCCSPL ( I ) ,RX, 13HCCSPL ( 1 + 1 ) /> 

WRITE  (6,6100  )(( I ,CCSPL ( I ) ,CCSPL(  1 + 1 )),  1 = 1 1 129,2  ) 
COMPOTE  THF  MONTE  CARLO  AVERGES 


<*00  IE  (NSAMP.LE.30)  GO  TO  10 
I =NSAMP— 30 
T = I 

XERRNL  = 08SSl)M(CERRNL(31  ; I ) ) / T 
XESONL  = OPSSOM ( CESONL ( 3 1 ; I ) ) /T 
XCSNL=08SSUM  (CCSMi.  ( 3 1 ; I ) ) / T 
XERRPL=OPSSOM ( CERRPL ( 3 1 : I ) ) /T 
X E SC PL =08 S SUM ( CESQPL ( 3 1 : I ) ) /T 
XCSPL=Q8SSUM (CCSPL ( 3 1 ; I ) ) / T 
X F VARNL =XESONL-XERRmLmxERRML 
XOBNL=DREPS 

IE  (XEVARNL.GT  .EPS(  1)  ) XD8NI.  = 10 .0*  A LOG  10  ( XE  V ARNL  ) 

XEVARPL=XESOPL-XERRPL*XFRRPL 

XUBPL=DBEPS 

IE  (XEVARPL.GT  .EPS(  1 ) ) XDBPL= 1 0 . # ALOG 1 0 ( XE V AR PL ) 
PRINT  THE  MONTE  CARLO  AVERAGES 


r 


% 

■ i 


I itf 


WRITE  (6,6050)  NMC , X ERRNL , XERRPL , XE SONL , XESQPL t XFV ARNL . 
1 XEVARPL  ,XDBNL  ,XDBPL  ,XCSNL  ,XCS PI- 
6050  FORMAT ( 1 H 1 , 2 5X , 3 0HM(  INTE  CARLO  SUMMARY  STATISTICS// 

1 32X , 1H ( , 14, 14H  SAMPLE  PATHS)// 

36X, 16HN0NLINEAR  F I L T ER , 5 X , 1 5HPH ASE  LOCK  LOOP// 

14X , 14H*##VARI ABLE*#* ,7X, 12H***VALUFS***,9X , 
17H***VALUES***//14X, 13HAVERAGE  ERROR,  9X, 
E15.8,6X,E15.8/ 


1 

2 

3 

6 

5 

6 

7 

8 
9 


i.  r ’ < ’ * * » ' ' * r r jl  -r  ? i jiih  v L.  ^ Fr\  f f 

E15.8,6X,E15.8/ 

10X , 21HAVERAGE  SQUARFD  ERROR , 5 X , E 1 5 . 8 , 6X , F 1 5 . « / 
1 4X , 1 4HERR OR  V AR I ANCF , 8 X , E 1 5 . 8 , 6X , E 1 5 . 8 / 

14X,14H  VARIANCE  IN  f)B  , 8 X , 6 1 5 . 8 , 6X  , E 1 5 . 8/ 
11X,19HAVERAGE  CYCLE  SL I PS , 6X , E 1 5 . 8 , 6X ,E 1 5 . 8 ) 


t 
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r 


best  mm  c:i 


4* 

0 


6100  F0HMAT(20X,I4,3XtEli>.6,4X,E15.8) 

GO  TO  10 
500  STOP 
FND 

SOHK  OUT  I NF  VPROP(A,I) 

REAL  A ( 16640  ) 

INTEGER  I 

INTEGER  MDl»MD2»M03»Mn4tMI)5  ,M|)6,Mn7» 

1 MD8,MD9,MD10,MD11 » MO  12, MO  13, MO  14 

COMMON  /CPROP/  M01  ,M02,M03,M04,M05,M06,MD7, 

1 MOR,MD9,MD10,M011 ,MD12,M013,M014 

IF  (I  .GT.O)  GO  TO  10 

a<md2;mdi ) = a < l ; MO l > 

A ( M04; MD3 ) =A ( 1 ;MD3 ) 

10  A(mq6;MD5)=A{ l;M05> 

A(md8:mo7)=A(1;md7) 

A(M010;MD9)=AI 1 ;MD9) 

A(MD12;M01 i )=A( l ;mdi l ) 

A(*D14:MD13)=A( 1 ;M0 13 ) 

R - TOR  N 
END 

SU8R  OUT  INF  M002PI  ( A ) 

REAL  A< 130) ,X( 130) ♦ Y ( 130) 

M I T rt T ( 130) 

COMMON  /PROB/  T WOP  I, PI 
M 1 : 1 30 ) =P  I 

10  "T(  1 : 1 30 ) = ( A ( 1:130)  ,GT.X(  1; 130)  ) 

IF  (O0SCNT  ( BT  ( 1 5 130)  ) .FO.O)  Gil  TO  20 
Yfl  : 130)=A ( 1 : 1 30 ) — T WOP  I 

*-(1:1  30  ) =05  vm  ASK  ( Y(  1 ; 130) , A(  1 : 130)  ,HT  ( 1 : 130)  : A ( 1 : 130)  ) 
GO  TO  10 
2 0 x ( 1 ; 130 ) =-PI 

3r  RT(  1 ; 130)  = ( A ( 1 : 130)  .LT.X(  1 : 130)  ) 

IF  (08SCNT(6T( 1:130) ). E0.0)  RETURN 
Y ( 1 ; 1 30  ) = A ( 1 ; 130  l+TWOPI 

Ail; 130) =08 VM ASK ( Y (1 ; 1 30 ) , A ( 1 : 1 30  ) , RT ( 1 ; 1 30  ) ; A ( 1 ; 1 3 o ) ) 

30  TO  30 
END 

REAL  FUNCTION  VAVG ( AV  » X , I ; # ) 

: F SCR  I P TOR  A V»  X » V AVG 
INTEGER  I 

= eai.  xi 

X I = I 

'/AVG  =(  ( X I - 1 . )*AV  + X)  /XI 

RE  TURN 
END 

SFAL  FUNCTION  RNF(INIT) 

INTEGER  INIT.K, KD.Ml t M2 tMl ,N2, NT, MP 
REAL  T 1 *T2 

COMMON  /RNUM/  K,K0,M1,M2,N1,N2,NT,MP,T1 tT2 

IE  ( INI T.EO.O)  GO  TO  10 

'<1  = 102943 

N2=  185617 

N 1 *2447  34 

N2= 1 58551 

Ml*?44734 


- K 
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irjii  -It  -*. 


M2=158551 
T1 *2 •+* ( -18 ) 

T2=2.**|-36) 

MP=2**IR 
10  K=M2*N2 
KI)=K/MP 
NT=K-KD*MP 
K = Ml«N2-*-M2*Nl  + KD 
KL)=K/MP 
N1=K-KD*MP 
N2  = NT 

RNF=N1*T1+N2*T2 

RETURN 

END 

REAL  FUNCTION  GAUSS ( INIT»S0»XM) 

INTEGER  INITfI*J 

REAL  SD,XM,TW0PI,X,XRlfXR2 

COMMON  /GNUM/  I,J,XR2 

COMMON  /PR08/  TWOPI 

IF  ( INIT.EO.O)  GO  TO  10 

1 = 1 

J=1 

10  IF  (J.NE.l)  GO  TO  20 
J = 2 

XR 1 =R NF ( I ) 

1=0 

XR2=RNF { 0 ) 

X = S0RT(ABS(-2.*AL0G(XR1  ) ) ) 

XR?=TW0P I#XR2 
XR1=X*SIN( XR2 ) 

XR2=X*=CnS I XR2 ) 

GAUSS=XR1*S0+XM 
RETURN 
?0  J=1 

GAUSS=XR2*S0+XM 

RETURN 

END 

SUBROUTINE  ST  AT E ( MC , SAMP , X 1 ,Z1,Z2  ) 

INTEGER  MC  * SAMP  t I N I T 

REAL  XI tZ 1 ,Z2,DEV1 , OE V2 T DE V3 ,0F V4 , X2 
COMMON  /STA/  0EV1,0EV2,DEV3,DEV4. INIT 

COMMON  /PROS/  TWOPI ,PI ,ALP1 10,0ELF ,022C,Y1EST ,Y2EST, All , A22, 
1 C0NST,DEL»FTC»PIDEL»P110«RDFL*RX,QQ,O22»M0*ND 
IF  (MC.NE.O)  GO  TO  10 
OE  Y 1 =SQRT (All  ) 

DEV2  = S0R  I"  ' A22  ) 

0EV3=SQRT (ROEL ) 

OEVA  = SORT (022  ) 

INI T = 1 
RETURN 

10  IF  (SAMP.GT.O)  GO  TO  20 
X l=GAUSS  < INIT,0EV1 tY]  ESV  ) 

I N I T =0 

X2=GAUSS (0,0EV2,Y2EST  ) 

RETURN 
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INTEGER  MC,SAMP,MD,ND,MND,MNDl 

INTEGER  MDl,M02,M03,M04,Mn5,MD6,Mn7,Mnft,MD9,MD10,MDll,Mt)l2,Mni'» 
MD9,MD9,MD10,MD11 , MO  12 , MO  1 3 , M014 
REAL  Z1»Z2»SHAT»CHAT,TNLF 

INTEGER  I,I1»J»J1»J2,K,NC,NT,NTERM,NTERM33 


i IF  (SAHP.LE.I)  GO  TO  30 
X 1 *X 1+X2*DEL 
X?=GAUSS(0,0EV4,X2) 

Z1=GAUSS(0,DEV3,C0S(X1) ) 

7?=GAUSS(0,DEV3,SIN<  XI ) ) 

RETURN 

END 

SUBROUTINE  PLL (MC,SAMP,Z1 , Z 2 , X 1H AT  , PF  1 1 ) 

INTEGER  MC , S AMP 

REAL  Z1 ,Z2,XIHAT,PF11 ,DEN , S INF  1 ,COSF  1 , 

1 X2HAT,PD11»PD12»PD22»PN11,PN12,PN22 

COMMON  /PLO/  0EN,X2HAT,PN11  ,PN12,PN22 

COMMON  /PROB/  T WOP  I , P I , AL P 1 1 0 , DE L F , 022C , Y I E ST , Y2E ST , A 1 1 , A22 , 
1 CONST f DEL » FTC, P IOEL , PI 10,RDEL,RX,00,022 
IF  (MC.LE.O)  RETURN 
IF  (SAMP.LE.O)  GO  TO  20 
IF  (SAMP.LE.l)  GO  TO  10 
X1HAT=X1HAT+DEL*X2HAT 

P011=(RDEL*(PN1 1+2.0*PN12*0EL  )-PN12*PN12 
1 *DEL*DEL)*OEN+PN22*OEL*DEL 
PD12=PN12*  (RDEL-PN12«0EL  ) *DEN+PN22*DFL 
PD22=-PN12*PN12*0EN+PN22+022 
PN1 1=PD1 1 
PN 1 2=PD1 2 
PN22=PD22 

DEN=1.0/( PN1 l+RDEL  ) 
i PF1 1=PNI 1*RDEL*0EN 
SINF1  = SIN<  X1HAT ) 

C0SF1=C0S ( X 1HAT ) 

X1HAT  = X1HAT+DEN*(-PNI  1*SINF1*Z  1 + PN1  1*C()SF1*Z2  ) 

X2HAT  = X2HAT  + OFN*<  -PM2*S  INF  1 *Z 1 + PN 1 2*C0SF 1 #7 2 ) 

RETURN 

I X1HAT=Y1EST 
X2HAT  = Y2EST 
PN 1 1 = A 1 1 
P N22=A22 
PN 1 2=0 • 

0FN=1.0/(PN11+R0EL ) 

RETURN 

ENO 

SUBROUTINE  NL F ( MC , S AMP , Z 1 , Z 2 , SH A T , CH A T , TNL F ) 


POINTMASS  FILTER  FUR  TWO  DIMENSIONAL  PHASF  ESTIMATION 
PROGRAMMED  BY  KENNETH  0.  SFNNE  JUNE  24,  1976 


non  ooo  oon  noo  ooo  non  oooooono 


REAL  All  ,A22,CL,CNORM,LUNST  ,CR ,P I , P IDEL .022  ,T  , 
l TT,  Y1FST.Y2E ST, TEMP, TEMPI, TEMP2 

REAL  SI  (64)  ,S2(64)  , SIGMA  (64)  ,PSI  (256)  ,A(256) 

RFAl  I N < 8 > 

INTEGFR  DJNS 
INTEGER  JNS ( 1 6640 ) 

REAL  COSY ( 16640) ,DELJ( 16640) , J0( 16640) ,JN( 16640) , JN1  ( 16640) , 

1 JNA(33280),SINY(16640),SN1( 16640) 

DESCRIPTOR  DSl ,DS2,DSN1 ,DSN2 ,0JN ,DJNS ,D JNA 1 ,DJNA2 ,0DEL J ,0 JN 1 , 
l D JNAENT  ,DJNABNT  ,()  JN  ANC  ,0  JN  A J 1,DJNAJ2,OSINY,OCOSY,OSN1A 

COMMON  /PROB/  TWOPI ,PI ,ALP1 1 0 ,DEL F ,022C , Y 1 E ST , Y2E ST , A 1 1 , A22 , 

1 C ONST, DEL, ETC, PIOEL, PI 10, RDEL,RX, 00,022, MD,ND 

COMMON  /NLFC/  NC,NT,NTERM,NTERM33,Sl,S2,SIGMA,PSI,A,CnSY,DELJ, 
1 JO, JN, JNS,SINY 

COMMON  /CPROP/  MO 1 , MO 2 ,MD3 , MD4 ,MD5 ,MD6 ,MD7 , 

1 MDB ,MD9,M010,M01 1 ,M0 1 2 , MO 1 3 ,M0 14 
IF  (SAMP.LE.O)  GO  TO  100 


SAMPLE  path  update  takfs  place  hfrf 


SET  CLOCK 

CALL  03CL0CKS ( T , TT  ) 

EVALUATE  SENSOR  T F R v S 

0SM  = Z1*DS1 

DSN2=Z2*DS2 

OS'  1 = DSN1+0SN2 

D S N 1 = V F X P ( 0 S N 1 :QSN1 ) 

SMI  M01  )=0. 

PROPAGATE  SENSOR  TERMS 
CALL  VPROP ( SN 1 , 0 ) 

SCRAMBLE  THE  JN  TO  OROFR  FOR  J(N+1) 


OJNA  1=08VGATHR ( 0 JN , 0 JN S J 0 JN A 1 ) 

FORM  THE  INTERPOLATED  MATRIX 

DJf'  =DJNA2-DJNA] 

OJN 1=0DELJ*0JN 
0 JN A 1 * DJNA 1+D JN1 

COPY  THE  END  COLUMNS 

DJNA  ENT  s»DJ  NAB  NT 

INITIALIZE  CONVOLUTION  (A(0)=1) 

0 JN 1 =Q JNANC 
C 

-y  .• 
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j 


c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


CONVOLUTION  LOOP 

J1=NC 
J?  = NC 

0(1  10  I = 1 1 NT  ERM 
.11  = J1  + MD1 
J2=J2-MD1 

ASSIGN  DJNAJ1 ,JNA( Jl;MN0) 

ASSIGN  DJNAJ2, JNA(J2;MN0) 
DJN=DJNA J 1+0 JNA J2 
D JN  = A ( I ) *DJN 
10  I)JN1  = DJN1+DJN 

MULTIPLY  BY  SENSOR  TERMS 

DJN1=DSN1A*DJN1 

GET  NORMALIZATION  CONSTANT 

CNQR  M = Q8S  SUM ( 0 JN  1 ) 

CN0RM=1 .O/CNORM 

TRANSFER  THE  NORMALIZED  DENSITY 
DJN=CN0RM*DJN1 
CUMULATE  ESTIMATES 
0JNAl=f)SINY*0JN 


* 


C 

C 

c 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 


ShaT=0BSSUM(DJNA1) 
I)JNA1=OCOSY«DJN 
CHAT=08SSUM(DJNA1 ) 

TIMEOUT 

CALL  D3CLUCKS(TNLE,TT  ) 
RETURN 


SAMPLE  PATH  INITIALIZATION  TAKES  PLACE  HERE 
10U  IE  (MC.LE.O)  GO  TO  200 

TRANSFER  THE  INITIAL  DENSITY  FOR  MEW  SAMPLE  PATH 


ooo  nno  nonnonooo 


GLOBAL  INITIALIZATIONS  OCCUR  HERE  FOR  THE  ENTIRE  RUN 


DETERMINE  THE  number  OF  CONVOLUTION  POINTS 

NTERM=(ND/2. )*SQRT ( SO . *022 ) /P IDFL+O . 5 

M01=MD+1 

M02=M01+1 

M0B=M01 *2 

MD4=MD3+1 

MD5=MD3*2 

MD6=M05+1 

MD7=M05*2 

M08  =MD7+ 1 

MD9=MD7*2 

MD10=MD9+1 

M01 1=M09*2 

MD12=MD1 1 + 1 

MD1 3 = MD1 1 *2 

M01A=MD13+I 

MND=MD1*ND 

MMO I =MND+ 1 

MT6RM33=MD1*NTERM 

NT=2#NTERM33 

NC=NTFRM33+1 

SET  UP  THE  VECTOR  DESCRIPTORS  FOR  THE  UPDATE  FUNCTIONS 

ASSIGN  DS1 ,S1 ( 1 JMO) 

ASSIGN  DS2,S2< I JMO) 

ASSIGN  DSN1 ♦ SMI ( 1 ;MD) 

ASSIGN  DSN2,JNA( I JMO) 

ASSIGN  DJN,JN(  UMNO) 

ASSIGN  DJNS, JNS( 1JMND) 

ASSIGN  DJNA1 , JNA< 1 ; MND ) 

ASSIGN  DJNA2, JNA ( 2JMND ) 

ASSIGN  DDEL J » DEL  J I 1 JMND ) 

ASSIGN  DJN1 , JN1 ( 1 JMNO) 

ASSIGN  DJNAENT , JNA ( MNDI JNT ) 

ASSIGN  DJNABNT, JNA( 1 JNT) 

ASSIGN  OJNANC t JNAINC JMND) 

ASSIGN  DSN1A, SNI ( I JMND) 

ASSIGN  DSINY,SINY< 1 JMND) 

ASSIGN  DCOSY,COSY( 1 JMND) 

PHASE  VARIABLES 

DO  210  1=1, MO 

SIGMA( I )=PI*( (2,*I-1.)/MD  -1.) 

COSY ( I ) =COS ( S I GMA ( I ) ) 

SI  NY ( I )=SIN(SIGMA( I ) ) 

SKI  )=COSY(I  ) /RDEL 
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r 


I 


V 

ft 


I- 


c 

c 

c 


c 

c 

c 


c 

r 

c 


c 

c 


c 


c 

c 

c 


210  S2< I)=SINY( I >/RDEL 
SINYIM01 )«0. 

CnSY(MDl)=0. 

CALI.  VPR0P(SINY,0) 

CALL  VPROP(COSY.O) 

PHASE  RATE  VARIABLES 

DO  220  1=1, ND 
220  PSI ( I )=PIDEL*< (2. *1-1. ) /NO  -1.) 

SETUP  THE  TRANSFER  MATRIX 

DO  24 0 J = 1 » NO 

J1=MDD( ND-l-NTERM+J ,ND ) *M01 
J2=( J-l ) *MD 1 
on  2 30  1 = 1. MO 

I 1 = Jl+l+MOO (MD+MD/2+32-< 135-NTFRM+J ) /4+I ,M0) 
I 2=J2+I 

230  JNS ( I 2 ) = 1 1 
J1=J2+MD1 

240  JNS ( J 1 ) =JNS ( J2+1 ) 

SETUP  INTERPOLATION  ‘ATRIX 

INI  1 )=0.fi75 
I N (?) =0. 625 
!M  3 ) =0 . 37  6 
IN| A) =0.125 
I I 5 ) = I N ( 1 ) 

I N ( M = I M ( 2 ) 

I " (7)  = I M ( 3 ) 

I N ( ? ) = I N ( 4 ) 

J = M 0 0 ( M T E R M , 4 ) 

DO  ?45  1=1,4 
I 1=  ( 1-1  )*Mf)l  + l 
,11  = 1+  4—  J 

T = I N ( J 1 ) 

245  ,’ELJ  ( I 1 :MD1  )=T 

CALL  VPROP ( DEL J » 1 ) 

EVALUATE  CONVOLUTION  TERMS  A(I) 

DO  ?«0  1=1 .NTFRM 
T=  I 
TT  = ND 
TEMP=T/TT 

temp=const*temp#temp 
A ( I ) =0. 

IF  ( TFMP.GT.-47)  A( I)=FXP(TFMP) 

7*0  CONTINUE 

CONSTRUCT  THE  A PRIORI  DENSITY 
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UPWMUHI,,!*; 


<•*»  iNmim 


M 


CN0RM=1.0/(TW0PI*SQRT( A11*A22) ) 

CL—0.5/A  22 
SI=-0. 5/All 
DO  290  1=1, MO 
I 1 = 1 

CR=SIGMA ( I I-Y1EST 
CR=CR*CR*S1 
00  290  J= 1 , NO 
TEMP=PSI ( JI-Y2EST 

J0(  1 1 )=EXP(TEMP*TEMP*CL+CR)  *CNORM 
290  1 1 = 1 1 + MD1 
C 

C WRITE  OUT  PARAMS  OF  NLF 

C 

WRITE  (6,6000)  MD,ND,M01,ND 

6000  FORMAT ( 1H0,26X,27HP0INT  MASS  NONLINEAR  F1LTER//1H  , 

1 28X, 2 4H VERSION  2,  COOED  6/27/76//1H  , 

2 18X, I 3, 1HX, I3,25H  DENSITIES  REPRESENTED  BY  ,I3,1HX,I3) 
WRITE  (6,6001)  NTERM,(A(I) ,I=1,NTERM) 

6001  FORMAT ( 1H  , 33X , 7HA ( 1 ) -A ( , I 2, 2H  ) / ( 1 X , 5E 1 5 . 8 ) ) 

RETURN 

C 

C 

C 

END 

0.0  0.0  -3.01  0.1  0.1  200  130  6^ 

0.0  0.0  -1.79  0.1  0.01  200  130  64 

END  OF  FILE 


C.  THREE-DIMENSIONAL  CDC-STAR-lOO  PROGRAM 


I 


* 


63 

64 


PROGRAM  MAIN(  INPUT  , HUT  PUT  , PUNCH , T APE  6=0UT  PUT  »TAPF5=INPUT  ) 
DIMENSION  XOAT  ( 130 , 1 0 ) , XHAT  ( 4 ) , Y 1 ( 35)  *Y2(  135)  * 

*EXPON(  35)  ,EXDON<  17 ) , Y3 ( 17 ) , PNF ( 3 , 3 ) , YA ( 35 > ,EXP33 < 16, 16 > , 

*0(3,3) , PBAR (3,3), PN (3, 3), AN (3) , F I 3 , 3 ) , PDUMMY (3,3), PDUM Y2 (3,3) 
COMMON/GN/ JGAUSS ,XZZZ ( 2 ) 

REAL  COSY ( 24576 ),SINY( 24576) ,SN 1(24576 ),SN2( 1536), JN( 24576 ) , 

* JN 1 (26112) ,JNA( 54272) ,nELJ( 26112) , SI ( 16) ,S2( 16) ,YB< 24576) ,D(2000) 
INTEGER  JNS( 1536 ) ,JNF ( 3072) 

BIT  B3I26112) 

OF  SCR  IP  TOR  DB,KJNF,DJNA,KJNS»njNl,DC 
ASSIGN  DB,JN(l;16> 

ASSIGN  KJNF,JNF( 153072) 

ASSIGN  DJNA,JNA( 1J49152) 

ASSIGN  KJNS»JNS( 1:1536) 

ASSIGN  0JN1,JN1( 1:26112) 

ASSIGN  DC , JNA (1;17) 

JG AUSS-0 

REA0(5,63)Y3EST,033C, ALF,GAM,NUM3 
FORMAT ( 4F 10,5,15) 

READ(5,64)Y1EST,Y26ST,  ALP1  10,DFLF,022C,NUM1  ,i\i(JM2  ,N02  ,NP3 
FORMAT ( 5F10.5,4I 5) 

P110=10.**( ALP110/10,  ) 

00=0220** ( .25) 

RX=(P110/( SORT ( 2.0)*00) )**(4,G/3.0) 

FTC  = S0RT(2.0)*RX**(  ,25)/00 
OELT=OELF*FTC 
022=022C*DELT 
R 1 1=RX /OELT 

P220=P110*S0RT (022C/RX) 

ALFD=ALF*DELT 
BET=1 .O-ALFO 

A 1 1 = 1 0. ** ( ( ALP1 10+GAM ) /10.  ) 

A22=P220 

P330=0.5*033C/ALF 
A33=?.0*P330 
033=033C*DELT 
0EVI  = S0RT( All  ) 

DFV2  = S0RT ( A22  ) 

DEV 3= SORT ( R 1 1 ) 

DE V4  = S0RT ( A33  ) 

DEV02=S0RT (022) 

DEV03=S0RT (033 ) 

KOUNT  = 1 
I S AMP= 1 
NS AMP  = 0 
SUMP  1=0.0 
SUMP2=0.0 
SUMP3=0.0 
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I 


y i 

k 

‘tJ 


OEL  SQ=DELT**2 
PI=A.*ATAN{1.  ) 

P 12=2. 0*P I 
PII)LT  = PI/OELT 
P I NV  = 1 . O/P I 
PI 2DLT=2.0*PI0LT 
III  = NIJM1 
U 2 = N IJ  M 2 
U3  = NIJM3 

IY2=U2/PI20LT*S0RT ( 50.0*022 )+. 5 

NTER  M=I Y2 

NTFRM1=NTERM+1 

NTERM16=NTERM*16 

NC=NTERMI6+1 

NT=2*NTERM16 

NTA1536=NT+1536 

R55=0. 5/R 1 1 

CL=-0.5/A22 

CM=-0.5/A33 

SI=-0. 5/All 

C 10  Y1,Y2  ANO  YA  a****#******:********************* 

EDG1=PI/U1 
EDG2=P 1 DLT/U2 
00  AO  1=1, NUM 1 
X = I - 1 
X=X/U1 

AO  Y1 (I ) =-PI +X*PI 2+E0G1 

00  50  I=1,NUM2 
X = I-1 
X=X/U2 

50  Y2 I I ) =-P IOLT+X*P 1 2DLT  +E0G2 
00  51  I =1 , I Y2 

X=  I 

X = X / 1 1 2 

51  Y 3 ( I ) =X 

00  5 5 1 = 1,  MUM  1 

COSY ( I )=C0S(Y1(  I ) ) 

5 5 S I MY ( I ) =S I N ( Y 1 ( I ) ) 

SKI;  16  ) =COSY  ( 1 ; 16)  /R1 1 
S2( 1 ; 1 6 ) =S I NY ( 1 ; 16  ) / R 1 1 
CALL  VPROP ( S I NY  , 1 ) 

CALL  VPROP (COSY, 1 ) 

XP=0.5*SQRT (A33  ) 

YY 3= (NUM 3-1 .01/2.0+1 .0 
I Y3=YY3 

00  60  1 = 1,  NIJM3 
X X = I -YY3 

60  YA ( I )=Y3EST+XP*XX 

C is**#######*****##*#*  OYMAMIC  EXPONENTIALS  ***$******************** 

I E ( I Y2. EO.O )G0  TO  153 
00  150  I = 1 , 1 Y2 

150  EX00N( I ) =EXP(-0. 5/022* (Y3( I ) **2 ) ) 

00  152  1=1, IY2 

EXP ON ( I )=EX00N( IY2  + 1-I  ) 

152  EXPON( I Y 2+1+1 ) =EXOON ( I ) 


x -mcmw  ' 


"‘r-’r*'1 1 nil 1 wru 


EXPONI IY2+1 >=0 
1YY=2*I Y2+1 
LTER*2=0 
nn  720  1=1,16 

00  720  K = 1 , 16 
\M)M=(K-I  )*XP+ALF0*(YA(  I ) 
\N0M  = -0.  5/Q33*XNUM**2 

EX  P 33 ( I ,K ) =0 . 0 

1 F (XNUM.LT.-27.  )G0  Til  720 
EX P 33 ( I , K ) =EXP(  XNIJM) 
CONTINUE 

LTERM=0 

lTERM16=LTERM*16 

LTERHl=LTERM+l 

LC=LTERM16+1 

NS=NTA1536*LTERM+NC 


00  339  K = 1 ,16 
DP  339  N= 1 , 16 

00  339  J =NT ERM 1 , I YY 

1 = 1 + 1 

0(1 )=EXPON(J)#FXP33(N,K) 

S it  # 4:  J*:}:***####:!:*#!*#  5*  £$!,>•.  INITIAL  DENSITY  * * * ♦ 3 

I JK  = 0 

00  160  K = 1 , 1 6 
ZZZ=CN*(YA(K)-Y3EST)**2 
00  160  J= 1 ,96 

YYY=ZZZ+CL*  ( Y2  ( J J-Y2FST  )**2 
CO  160  1=1,16 
! J*  = I JK+ 1 

XXX=YYY+SI*(Y1 ( I )-Y16ST)**2 
IFIXXX.LT.-27. )G0  TO  159 
JM(  I JK ) =EXP ( XXX  ) 

SO  TO  160 
JM  I JK  ) =0 .0 
CONTINUE 

a##*##  ##!*#*>*>*  INTEGER  ARRANGEMENT  *$*#*=**4:#*!* 
DO  225  1=1,26112 
“31  I ) =rt • 1 ' 


00  300  K=  1 , 1 6 

00  300  J = 1 ,96 

1 = 1 + 17 

3 3 ( I ) =B • 0 • 

00  320  K=1 ,3071 ,2 
J\F ( K ) = ( K-l )*« 

JN E ( K + 1 ) = JNF ( K > 


J1=0 

00  332  K = 1 , 1 6 

00  332  J= 1 ,96 
1=1  + 1 

I l*Jl+MOD(23-( J-l ) /6 , 1 6 ) 
I ) = Il 
J1=J 1+32 

DELJ( 1 5 17 )=11./12. 
0ELJ(18;17)=0.75 


DELJ<35;17)=7./12. 

DEL  J ( 52;  ?. 7 ) =5.  / 12. 

DELJ<69;17)=0.25 
HFLJ<«6; 17)=1./12. 

CALL  VPROPl(DFLJ) 

11  CONTINUE 
Kl  IUNT  = 1 

CALL  GAUSS ( JSEED , DEV  1 , Y IE  ST , X 1 ) 

X DA  T ( K DIJNT  , 1 )=X1 
CALL  GAUSS ( JS FED , DE V 2 , Y2E S T , X 2 ) 

CALL  G AUSS ( J SEED »0E V4 , Y3EST  ,X3) 

AC0S=X3*C0S (XI ) 

AS I N = X 3*S I N ( X 1 1 

CALL  GAUSS( JSFED,DEV3, ACOS,Zl ) 

CALL  GAUSS( J SEED , DE V3 , AS  I N , Z2 ) 

GO  TO  470 
450  CONTINUF 

X1=X1+X2*DELT 
X D A T ( K OUNT  , 1 ) =X 1 
CALL  GAUSS ( JSEF0»DEVQ2*X2,X2 ) 

X3  = X3*BET  + AI_FD 

CALL  GAUSS( JS EED , DE VQ3 , X3 , X3 ) 

AC0S=X3*C0S(X1 ) 

AS  I M = X 3*S I N ( X 1 ) 

CALL  GAUSS! JSEED, 0EV3,ACPS,Z1  ) 

CALL  GAUSS! JSEED, DEV3,ASIN,Z2) 

XP=0.5*SQRT (A LOSS) 

DO  600  1=1, NUM 3 
X X = I -YY3 

YA  ( I ) =X3EST+XX*XP 
600  CONTINUE 

DO  730  1=1,16 
DO  730  <=1,16 

X NtIMs  ( K—  I ) ~XP  + AL  FO*  ( Y A ( I ) — 1 . ) 

XIMUM  = -0,5/Q33#XMUv'#*2 
EXP33 ( I , K )=0.0 
IF(XNUM.LT.-27.)G0  TO  730 
EX  P 33  ( I ,K)=EXP(XIMIJM) 

730  CONTINUE 
1=0 

DO  3A0  K= 1 , 1 6 
DO  340  N = 1,16 
DO  340  J = NT  ERM 1,1 YY 
1 = 1 + 1 

340  D( I )=EXPON( J)*EXP33(N,K) 

470  CONTINUE 

i C #####$#**#*#*«*>*#***!)!#«#  SENSOR  FUNCTION  *#«#***#«#**$***#**>:<##**  = 

•4  CALL  Q3CL0CKS ( T ,TT  1 


S1(1:16)=Z1*S1( 1 ; 16 ) 
S2(l?16)=Z2*S2(l;16) 

SI ( 1 ; 16 ) =S1 ( 1 ; 16 ) +S2 ( 1:16) 


J 


I 


[ ‘ i 

| H ■ 

I 

f t? 
» 

i \ ■ 

1 1; 
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CALL  VPROP ( SN2  *0 ) 

SN1 (J;1536)=SN2( 1 *1536) 

S2 ( 1:16)=-R55*YA(K)*YA(K) 

SN2<  1 ; 16  >=VEXP<  S2 ( 1 5 16) ; SN2( 1 J 16  ) ) 

CALL  VPROP ( SN2 ,0 ) 

SN1 ( J; 1536 )=SNl( J: 1536)*SN2< 1; 1536) 

500  J = J+ 1 5 36 

C **$$**#$##***$*#*£*  MAIN  LOOP  STARTS  *#*****##*$*$#$#*#**«*#♦***** 

CALL  08 VX T OV ( X 1 02 ’ , 0 , K JNF , 0 , DH , 0 , 0 JN A ) 

CALL  08VXT0V(X  *02'  ,0,KJNS,  0, Of. , 0 , 0 JN  1 ) 

JNA (1:26112) =JN1(2; 26 112) -JN1(1;26112) 

JNA  ( 1:26112  )=I)EL  J ( 1:26112)*  JNA(  1:26112) 

JN1 ( 1:261 12 )=JN1 ( 1 ;261 12)+JNA< 1 ;261 12) 

JNA ( 1 ; 24  576 )=08VCMPRS( JN1 ( 1:26112) , R 3 ( 1:26112) : J N A ( 1 ; 245  76 ) ) 

JN 1 ( 1 ; 24576 )= JNA ( 1:24576) 

J = 1 
1 = 1 

DO  510  K = 1 * 1 6 
N= I +NTERM16 

JNA(N; 1536) =JN1( J: 1536) 

JNA ( I ; NTERM 1 6 ) = JN A ( 1 + 1536: N TERM  16) 

JNA ( N+l 536 ; NTERM 16 ) =JM 1 ( J;NTFRM 16) 

J= J+ 1 536 

510  1 = 1 +NT  A 1 536 

N=1 
11=0 
JK=1 

00  700  1=1,16 
JMIN!  1536  ) =0.0 

00  690  K= 1 , 1 6 
J 1 =NS  + 1 1 
J2=NS+I 1 

L>0  680  J = l, NTERM  1 

JM( 1 ; 1536 ) =JNA ( J1 ; 1536) +JMA ( J2: 1536) 

JN(l: 1536)=JN( 1 : 1536 )*0( JK ) 

JN 1 ( M;  1536 ) = J N 1 ( N:  1 536) +JN ( 1 ; 1536) 

JK  = JK+ 1 
J1=J1+16 
680  J2=J2-16 

690  CONTINUE 
N=N+ 1536 

700  I 1 = 1 1+NTA1536 

JN1 ( 1 : 24576 )=JN1 ( 1 ; 24576 )#SN1  ( 1 : 24576) 

CN08M=SUML0G ( JN1 ) 

IE(CN0RM.LT.1.0E-20)CN0RM=1.0 
CNORM=l ,/CNORM 

J N ( 1 ; 2 45  76 ) =CN0RM*JM1 ( 1:24576) 

JNA ( 1 : 24576 )=COSY( 1 :2457)*JN( 1:24576) 
r.HAT  = SUMLOG ( JNA ) 

JNA ( 1 ; 24576 )=SINY( 1 ; 24576) *JN( 1:24576) 

SHAT=SUML0G( JNA) 

CXH  A T = A T AN2 ( SHAT , CHAT ) 

J=1 

DO  343  K = 1 , 1 6 
YR(  .1  : ) 536  ) =YA  f K » 
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ii 
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343 


J= J+1536 

JNA ( 1 ; 24576 ) = YB  ( 1 ; 24576 )*JN( 1 ; 24576) 

X3EST=SUHL0G( JNA) 

JNA  ( 1 ; 24576 )= JNA ( 1 ; 24676 )*YR( 1 ; 24676) 

ALOSS  = SUMLOG (JNA  ) 

AI.OSS  = ALOSS-X3EST*X3HST 
c ifc**#*###*!!*************#***  MAIN  LOOP  ENOS  #4#$ YS#**  *************** 

CALL  03CL0CKS(TNLF,TT ) 

WRI TE (6,201 )KOUNT ,X1 ,X2,X3,Z1 , Z2 , C XH AT , X3E ST  , A LOSS 
201  FORMA  T ( I 5 , 1 X , 1 P3E 1 4 . 6, 4X , 1 P2E 1 4. 6, 4X , 1 P3E 1 4 . 6 ) 

WRITE(6,8880)TNLF 
8880  FORMAT ( 1PE12.6) 

I F ( KOUNT ,EQ.N02)G0  TO  505 
KOUNT  =K0UNT+1 
GO  TO  450 
505  CONTINUF 
SUMP  = 0 , 0 
SUMC=0.0 

X OAT (KOUNT, 2 )=CXHAT 
XDAT ( KOUNT , 3 ) =ALOSS 
XOAT (KOUNT ,4)=X3EST 
DO  1501  1=31, N02 

XD=ABS (XDAT( I , 1 )-XOAT  ( I ,2)  ) 

1498  CONTINUE 

IF(XD.GT.PI  ) GO  TO  1499 
GO  TO  1500 
1699  X D = X 0— P I 2 
GO  TO  1498 

1500  SUMP  = SUMP  + Xf)**2 

1501  CONTINUE 
H=N02— 30 
sump=sump/h 

XNSAMP=NS AMP 
XAA  = XNSAMP+1 .0 

SUMP1  = ( SUMP+XNSAMP*SIJMP1  ) / X A A 
D SUMP  1 =AL0G10 ( SUMP  1 )vl0.0 

00  1601  1=31, N02 

XD  = ABS ( XDAT ( I , 1 ) -XOAT ( I ,4  ) ) 

1698  CONTINUE 

I F ( XD.GT .P I ) GO  TO  1699 
GO  TO  1700 

1699  XD  = XD-P I 2 
GO  TO  1698 

1700  SUMC=SUMC+X0**2 
1601  CONTINUE 

SUMC=SUMC/H 

SUMP 2= ( SUMC+XNSAMP«SUMP2 ) /xaa 
0SUMP2  = AL0G10( SUM P2  1*10.0 

WRITE  (6, 1511  ) N03  , SUMP  1 , OSIJMP  1 , SUMP2,OSUMP2 
1511  FORMAT* I10,1P4E14.6) 

MSAMP=NSAMP+1 

I F ( I SAMP . EO . N03 ) GO  TO  2200 

1 SAMP  = I SAMP+l 
GO  TO  11 


2200  CONTINUE 


STOP 

ENO 

FUNCTION  SUMLOG(A) 

REAL  A(26112> ,C< 12288) 

C( 1:12288 )*A<1 512288 )+A( 12289:12288) 
C( 1 56144) =C( 1;6144)+C<6145;6144) 

C< 1; 3072) =C< 1:3072) +C( 3073 5 3072) 

C( 1 5 1536) =C< 1 ; 1536) +C( 1537; 1536) 

C ( 1 5 768  ) =C ( 1 5 768 ) +C ( 769  5 768 ) 

C ( 1 ; 384) =c< 1 ;384)+C (385:384  ) 
C(1;192)=C(1;192)+C( 193:192) 

C( 1;96)=C( l;96)+c (97;96) 
C(1;48)=C(1;48)+C (49:48) 

C( i:24)=C( i;24)+c< 25:24) 
C(l;12)=C(l;12)+C(13;12) 
C(l;6)=C(l;6)+C(7;6) 

C ( 1 ; 3 )=C  ( 1 ;3  )+C  (A  ;3 ) 

SUMLOG=C ( 1 ) +C  ( 2 ) +C  ( 3 ) 

RETURN 

END 

SUBROUTINE  VPROP ( A, I ) 

REAL  AI2A576) 

IE( I .E0.2  >G0  TO  10 
A( 17; 16)=A ( 1 ; 16) 

A ( 33 ; 32 ) =A ( 1:32) 

A ( 65 ; 32 ) =A ( l;32) 

10  A ( 97 ; 96 ) =A ( 1 ; 96 ) 

A ( 193; 192 ) =A ( 1 : 192  ) 

A ( 38  5 5 384 ) = A ( 1 S 384  ) 

A( 769: 768 ) =A ( 1 ; 768 ) 

I F ( I . EQ. 0 ) R ETURN 
A ( 1537; 1536) =A ( 1 ; 1536) 

A( 3073:3072 ) = A ( 1;3072) 
A(6145;6144)=A( 1 ;61AM 
A ( 12289;  12288 )=A<  1 : 1 2288  ) 

RETURN 

END 

SUBROUTINE  VPROPl(A) 

REAL  A ( 26 1 1 2 ) 

A ( 103:102 ) =A ( 1 ; 102  ) 

A(205;204)=A{ 1 ;204) 

A ( 409:408 ) = A ( 1:408  ) 
A(817;816)=A(1;816) 

A ( 1633;  1632) =A( 1 ; 1632) 
A(2265;2264)=A( 1;2264) 
A(3265;3264)=A( 1;3264) 
A(6529;6528)=A( 1 ;6528) 

A( 13057; 13056 )=A( 1:13056) 

RETURN 

ENO 

FUNCTION  RNNF ( NS  »MODE ) 
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DIMENSION  NS ( 2 ) , NC (2) 

COMMON  /RN/  Nl,  M2,  MP,  Tl,  T2 
DATA  MI,  M2/244734,  158551/ 

C MODE=0  TO  CONTINUE,  OTHERWISE  RESTART  WITH 

C INTEGER  NUMBER  NS ( 1 ) *2**1 R+NS { 2 ) 

IF  (MODE)  10,  100,  10 
10  N 1 =NS ( 1 ) 

N2  = NS ( 2 ) 

T 1 =2 . ** ( -18  ) 

T2  = 2. **(-36) 

MP=2**18 

100  DO  200  1=1,2 

GO  TO  ( 110,120), I 
110  K=M2*N2 
GO  TO  190 

120  K=M1*N2+M2*N1+KD 
190  KD=K/MP 
200  NC ( I )=K— KD*MP 
N1=NC (2 ) 

N2=NC( 1 ) 

XN1 =N 1 
XN2=N2 

RNNF=XN1*T1+XN2*T2 

RETURN 

END 

SUBROUTINE  GAUSS ( JS  , SD , XM , X ) 

DIMENSION  NST ( 2 ) 

COMMON  /RN/  Nl,  N2 , MC , Tl,  T2 
COMMON  /GN/  J , XR ( 2 ) 

IF  (J)  10,  10,  20 
10  J = 2 

TW0PI=8.*ATAN( 1. ) 

NST ( 1 >=244734 
NST ( 2 ) =158551 
NST< 1 >=102943 
NST(2)=185617 
XR ( 1 ) =RNNF ( NST  , 1 ) 

GO  TO  35 

20  GO  TO  (30,40)  , J 
30  J = 2 

XR ( 1 ) =RNNF ( NST , 0 ) 

35  XR(2)=RNNF(NST,0) 

X1  = S0RT(ABS(— 2 . * ALOG ( XR ( 1 ) ) ) ) 

XR(2)=TWOPI*XR(2) 

XR ( 1)=X1*SIN(XR(2)  ) 

XR(2)=X1*C0S(XR(2)  ) 

X=XR ( 1 ) *SD+XM 
RETURN 
40  J=1 

X=XR ( 2 ) *SD+XM 

RETURN 

ENO 

1,0  0.01  1.0  1.4  lfe 

0.0  0.0  -3.0  0.1  0.01  16  96  130 

END  OF  FILE 


i 


) 
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